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ABSTRACT 


This report evaluates fortran coded out-of-core equation solvers 
thit solve using direct methods symmetric banded systems of simultaneous 
algebraic equations. The types of solvers studied were banded, frontal 
and column(skyline) solvers. Also considered were solvers that could 
"partition" the working area and thus could fit into any available core. 
Comparison timings are presented for several typical two-dimensional and 
three-dimensional continuum type grids of elements with and without mid- 
side nodes. Extensive conclusions are also given. 
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CHAPTER 1 INTRODUCTION 


The purpose of this report is to evaluate fortran coded out-of-core 
equation solvers that solve using direct methods symmetric banded systems 
of simultaneous equations such as the equilibrium equations generated by 
static finite element or finite difference discretizations of two- and 
three-dimensional stress analysis problems. This research was carried out 
under contract to the N.A.S.A. Marshall Space Flight Center for the purpose 
of improving the operation.:) efficiency of the TEXGAP computer program [1]. 
This report is directed to those individuals who develop, modify or rou- 
tinely use finite element computer codes. No attempt is made herein to 
review finite element methodology and the reader is assumed to be familiar 
with fortran programming techniques and the techniques for solving systems 
of simultaneous equations on large digital computers. 

There are 3 methods commonly used for solving systems of simultaneous 
equations on digital computers; direct methods, iterative methods and 
hybrid methods. Direct methods involve the reduction of the equations to 
an upper triangular form (this step is called forward substitution, reduc- 
tion or triangularization) from which solution is effected by simple back- 
substitution. Iterative methods require an approximate starting solution 
from which successively better approximations are obtained from algebraic 
equations using one of many procedures such as Gauss-Seidel [2] or Succes- 
sive Over-Relaxation [3]. Hybrid methods involve the combined use of direct 
and iterative methods. 

The earliest finite element programs developed in the late 50' s and 
early 60' s used the Gauss-Seidel iteration method to solve for the nodal 
point displacements [4]. Iteration methods were probably selected by these 
early developers because they were structural analysts familiar with a 
similar iteration technique called "Moment Distribution". These early 
iteration schemes quickly gave way to the direct method of Gaussian elimin- 
ation. Direct methods proved to be just as convenient to program and they 
were easier to use and gave accurate and reliable results for a large class 
of two-dimensional (2D) static linear elastic stress analysis problems. 
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At first, the direct equation solvers were programmed so that all the ^ 
coefficients fit within the available core storage. Earlier research also 
funded by NASA studied these solvers [5]. At the time of this earlier 
research (the mid 60' s), the available central memory core storage was 
generally either 32,000 or 65,000 words. Thus only 2-3,000 unknowns or 
degrees-of-freedom(DOF) could be handled using even the most efficient banded, 
symmetric solvers. As finite element methodology advanced, elements with 
more nodes and more DOF and three-dimensional (3D) applications were developed 
that required millions of words of storage, much more core than Is available 
on any computer system. This lead to the development of solvers that parti- 
tioned the equations into smaller pieces that would fit Into central memory 
core storage while the eliminations were being performed, and then these 
smaller pieces were written onto low speed tape, disk or drum storage to 
make room for the next set of coefficients. 

While the techniques for programming Gaussian elimination in-core 
are straightforward and fairly well standardized, those used for out-of-core 
solvers are highly variable. Another complicating factor Is that virtually 
every computer Installation has a different method of charging for the trans- 
fer of data from central memory to low speed storage. Some make no charge 
and others may charge more for these input/output (I/O) transfers than for 
the actual central processor(CP) calculation time. For example, the Uni- 
versity of Texas charges $0,004 per each 64 words (called a Physical Record 
Unit) but converts this charge to equivalent CDC 6600 time at the rate of 
$230 per hour. Thus, the charge rate amounts to 0.98 x 10*^ CP sec per 
word transferred. The Importance of this charge can be seen In the following 
example. A system of n*10^ equations at a bandwidth of b=10^ requires 
approximately nb^/2 » 5 x 10® operations, I.e., the total number of 
multiplications and additions to triangularize the system. This will require 
approximately 2 x 10^ CP seconds on a CDC 6600 assuming the average rate 
of computations Is about .25 x 10^ operations per CP second. Provided 
that the entire bandwidth block fits In-core, there will be 2nb = 20 x 10® 

I/O transfers (1 write and 1 read of each coefficient) or an equivalent 
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4 4 

charge of 2 x 10 CP seconds. Thus, the total charge will be 4 x 10 

seconds, fully one-half being charged for I/O transfers. 

The net effect of this is to make meaningful evaluation of oi’*-of- 
core solvers very difficult because the concept of operational eff iency 
means performing a given set of computations for the minimum computational 
expense (i.e., the actual charge) and not the minimum central processor 
time. Further complications which are not considered herein also arise when 
consideration must be given to obtaining reasonable priority to achieve good 
turn-around. 

In the next chapter is given background information on fortran coding 
techniques to solve sys*^ems of banded symmetric positive definite equations 
and specific attentior: is given to out-of-core bandad, frontal and column 
solvers. Chapter 3 briefly describes soriie important characterics of the 
particular solvers used in the present study. Chapter 4 describes and 
presents the results of the various numerical experiments carried out to 
evaluate these solvers. Chapter 5 gives rather extensive conclusions and 
recommendations for future research. 
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Background/In-Core Solvpr^ 


solvers w,1 now be made to establish the tenn1nolog.v used in this report 
In this section only in-core solvers are considered. 

Let the simultaneous equations be representated as 


n 

r 

j=i 




for f = 1,2,. 


( 2 . 1 ) 


where a^^ are the stiffness coefficients, b, the nodal point forces 
and the unknowns (DOF), and n is the total number of DOF. The standard 
el, m, nation procedure is to solve for x, in tenns of x x "V 
from the first equation. It is important to use the f?^st%ua«o"n to 

Solving the first equation for gives ^ 


This equation is now substituted into equations 2 thru n 
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and collecting 
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= 2,3,. 


evidence that^indicatewlord^ recognized body of 

error with systems of linear equations or reduces roundoff 

mixed finite element models. either displacement or 


Thus, the order of the system of equations is reduced from size n to n-1 
After n-1 such eliminations (Xj from the second equation, etc ), the 
equations are reduced to an upper triangular form that permits solution for 
the unknowns by backsubstitution. This direct Gaussian elimination is 
straight forward to code in the Fortran language as is shown in Table 2 1 


DO 300 K=1,N-1 
DO 200 I=K+1.N 

B(I)=B(I)-B(K)*A(I.K)/A(K,K) 

DO 100 J=Ki-l.N 

100 A(I.J)=A(I,J)-A(I.K)*A{K,J)/A(K.K) 

200 CONTINUE 
300 CONTINUE 
C 

C BACKSUBSTITUTION 
C 

X(N)=B(N)/A(N,N) 

I=N 

400 l-I-l 
SUM=0.0 

DO 500 J=I+1,N 
500 SUM=SUM-A(I,J)*X(J) 

X{I)=(B(I)-SUM)/A(I,I) 

IF (I.GT.l) GO TO 400 

Table 2.1 Symbolic Fortran Routine for Direct Gaussian Elimination 


n fh K either the symmetry 

or the banded nature of the equilibrium equations generated by finite element 

^thods As Illustrated in Figure 2.1 only those coefficients above the main 

diagonal and within the band need be stored because the Fortran coding given 

in Table 2.1 will not cause any nonzero terns to appear outside the band 
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during the elimination process*. Thus, it is necessary to store only n*b 
(b=bandwidth) coefficients instead of the x? required for a general system. 

To take advantage of this storage reduction, it is necessary to store the 
equations in a tri-diagonal form in which column 1 of the tri-diagonal matrix 
gives the main diagonal of the actual system of equatii ... column 2 the first 
off-diagonal term and column b the edge of the band in any row. The Fortran 
coding to effect reduction of the equations in tri-diagonal form is given in 
Table 2.2. The coding in Table 2.2 is essentially the same as that in 
Table 2.1, however, advantage is taken of the banding and symmetry to reduce 
the number of passes through the innermost loop (300) from n^3 in the gen- 
eral case to nbV2. 


DO 300 K=1,N-1 
LIM=MIN(K+MBAND,N) 

DO 200 I=K+1,LIM 
B(I)=B(I)-A(K,I-K)*B{K)/A(K,1) 

DO 100 J=I,LIM 

100 A(I,J-I+1)=A(I,J-I+1)-A(K,I-K)*A{K,J-K)/A(K,1) 

200 CONTINUE 
300 CONTINUE 

Table 2,2 Symbolic Fortran Routine for Tri-diagonal ized Equations 

There are many obvious Fortran improvements that can be made to the 
coding in Tables 2.1 and 2.2**, but these are not particularly relevant to 
the present study since the limiting feature of this coding is that it 
requires all coefficients to reside in-core during the solution. If we 
assume that the largest an array can be dimensioned is 4 x 10^ decimal 
locations, this limits the size of the problem that can be solved to say 


♦Reordering might cause the band to increase. 

**For example, the “DO” loops should be replaced by variable counters 
U-I+1, etc.), the variable subscripting should be with single subscripts 
and the innermost sum accomplished with temporaries. 
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n=800 and b=50. This method of solution Is then capable of handling one- 
dimensional (ID) continuum problems, 2D and some 3D truss and frame struc- 
tures and some small 2D continuum problems (a plane strain grid with 400 
nodes and a maximum bandwidth of 25 nodes). Clearly, In-core solvers 
have a very limited spectrum of applications. 


2.2 Out-of-Core Solvers 

It Is very difficult to characterize out-of-core solvers since a 
very Important feature is how the coeffic.ents are stored and how many 
I/O transfers are required. Herein is briefly described three Important 
types of out-of-core solvers; band, frontal and colurm solvers. For band 
and frontal solvers there is an Important dichotomy between those that have 
sufficient storage available to handle the maximum band or front and those 
that lack sufficient storage and thus must subdivide the band or front. The 
subdivision of the band or front will be called "partitioning". This Is not 
an entirely satisfactory name since partitioning is common used in matrix 
algebra to denote symbolic groupings, but It accurately conveys the Intended 
meaning In the present study. Out-of-core column solvers are intrinsicly 
partitioned by the manner In which "blocking" is done. These factors will 
be discussed more fully In the next 3 sections. 

2.2.1 Bandsolvers 

All out-of-core bandsolvers use logic essentially the same as that 
given In Table 2.2 for effecting the elimination process, but they differ 
significantly in how they perform the I/O. Figures 2.2 and 2.3 Illustrate 
the two extremes In terms of required core storage for bandsolvers that do 
not partition the band. 

^ The sliding block scheme illustrated In Figure 2.2 requires storage 
for b /2 coefficients and as a row Is eliminated it Is written to low 
speed storage and the equations are shifted up one position [6]. There are 
several drawbacks to this method: 
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(1) I/O is required row-wise and therefore transfers small records, 

(2) shifting takes place after each elimination which is a very 
wasteful procedure, 

(3) the procedure for assembly and also that for bringing new 
columns into core are very complicated, and 

(4) assembly and elimination cannot be performed simultaneously. 

The bandwidth block scheme illustrated in Figure 2.3 avoids all of 

the above inefficiencies and complexities at the expense of storage. A 
2 

total of 2b words are necessary because the first "b" unknowns are elim- 
inated before any I/O is done. This requires that the triangular portion 
in block II be present during elimination of block I. Actually, only 1/2 
of block II is needed but it is not generally convenient to take advantage 
of this space* so b^/2 of storage is wasted. A significant advantage to 
the bandwidth block procedure is that it permits (in fact encourages) the 
combination of the assembly of the global stiffness from the element stiff- 
ness and the forward elimination process. This is significant because 
it eliminates wasted I/O and limits the total transfers to n*b words (1 word 
for each coefficient). 

There are of course many other methods, but their differences are 
minor for purposes of the present study. 

Without partitioning it is necessary to fit some multiple (2) of 
b^ into ce ntral memory, thus if only 32,000 words are available then 
^^^J20M ="126. However, on CDC 7600s and many other large computers there 
is often 250,000 words available. In this case b£353, a very considerable 
number. 

For large bandwidths it is necessary to partition the band as illu- 
strated in Figure 2,4 for 3 partitions of the b^ block I, 1^ I^ & I^ 
respectively. This procedure implies core residence of 2 of the b^/3 
blocks at any time. Blocks I^ & initially are brought into core and 


♦Simply because the.coding becomes too complicated. 
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I.| is reduced through the end of I 2 , then partially reduced block I 2 is 
written out and is brought in and the elimination of is made on it. 

Ig is written and II^ is treated the same as is now complete 

and written and partially reduced and are brought back in and so on. 
Table 2.3 gives the total number of reads and writes (1 for each) for only 
the reduction. 

total number of reads & writes 

3 
5 
7 
9 
9 
9 
9 
7 

58 

Table 2.3 Block I/Os in Reduction 

Thus for this simple example, it is found that 58 I/Os of block size b^/3 
were required. Without partitioning there would have been 3 outputs of 
size b . The partioning scheme thus required 58b /3 t 3b^^ = 58/9 - 6 times 
as many I/O word transfers and 58/3 = 19 times as many I/O accesses. Note 
that block I must be written during the assembly process and then read and 
written once during the reduction process. Assembly and reduction cannot be 
combined in a bandsolver that requires partitioning. Obviously, partitioning 
should be avoided whenever possible. 

2.2.2 Frontal Solvers 

Frontal solvers are fundamentally quite different than either band 
or column solvers (which are similar) in that the nodes or DOF are assigned 
location numbers in the active front based on their appearance in the element 


block 


'2 

■3 

”1 

"2 

"3 

III, 

III, 
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being assembled. Thus, the ordering of the element stiffness matrices passed x 
to a frontal solver determines the size of the front. The nodal numbers are 
used only as reference numbers or labels with which to compare with other 
labels (e.g., Al, BIO, etc. are just as useful as 1, 210, etc.). 


To illustrate how a frontal procedure works consider the grid shown 
in Figure 2.5. As element 1 is assembled, nodes 1,2,3,8,9,12,13 & 14 are 
given (column or front) locations 1 .2,3,4, 5, 6,7 & 8. respectively in the 
"active" area of the band which is called the front. Since nodes 1,2 & 8 
are connected only to element 1 they are now eliminated from the active area 
in the usual manner and the reduced coefficients passed to an internal buffer 
area that is dumped to low speed storage when full. This leaves the active 
area with a substructure stiffness consisting of the line of nodes 12-13- 
14-9-3. Next element 2 is assembled with n^ nodes 4,5,10,15 & 16. The 


new nodes are assembled into the first available positions in the front; 
1.2.4,9 & 10, respectively. Nodes 3.9 & 14 are currently located in the 
front at locations 3.5 & 8. respectively. Since nodes 3,4 & 9 are complete 
they are now eliminated leaving the active front the nodal line 12-13-14- 
15-16-10-5. Table 2.4 summarizes the frontal locations and maximum front 
widths for all 6 elements. Zeroes (0) indicate that the node has not yet 
appeared, -1 indicates it has been previously eliminated, and 4* means 
that the node is in position 4 and will be eliminated because the current 
element is the last appearance of that node. Note that after an initial 
transient the front remains steady for regular grids. Note also that the 
frontal method inherently combines assembly and reduction. 


Frontal solvers are typically organized into 3 main sections; 
prefront, front and backsubstitution. The prefront constructs the table 
of first, intermediate and last appearances illustrated in Table 2.4. The 
front performs the assembly and reduction, and the backsubstitution the 
usual solution computation. All 3 sections are generally coded with 2 main 
nested loops. The outer loop on the elements and the inner loop on the 
number of DOF for that element. Of course other loops must be nested inside 
these two mam loops to perform the reduction and backsubstitution. It is 
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Figure 2.5 Example Grid of Q8 Elements 
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possible to code pre-prefronts that order the elements in an optional or 
sub-optional manner so as to minimize the front. This is similar to band 
minimizers and was not considered in the present study. Front and band 
minimizers are not useful for regular grids. 

The frontal method is especially efficient in handling grids with 
midside nodes as shown in Figure 2.5 and Table 2.4 where the maximum front 
is 12 nodes. The maximum bandwidth would be 14 nodes. In general if there 
are N elements across the width of a uniform grid of Q8 elements with 2 DOF 
per node, the front f=2(2N+6) but the band b=2(3N+5). Thus, in the limit 
as N-^ f/b^2/3. The frontal method is also much more efficient in 
storage and computations. The frontal method requires storage for K^=f^/2 
coefficients whereas the band method needs thus in the limit as 

K^/K|j-*^l/9, almost an order of magnitude. The number of computations 
is proportional to f^ and b^, respectively, thus the ratio of front to 
band approaches 4/9 as N-x». 

The frontal method loses its computational advantage for grids 
without midside nodes since f=b. However, it still retains a significant 
storage advantage of 1/4. 

2.2.3 Column or Skyline Solvers 

Column solvers are currently very popular among the coterie of 
developers of general purpose finite element codes (e.g., SAP [7,8J and 
ADINA [9]) and there have been many recent papers in the literature 
describing them [10,11 ,12, 13J. Column solvers are based on the skyline 
storage scheme illustrated in Figures 2. 6-2. 8 in which a nonuniform 
skyline is drawn covering the minimum number of nonzero coefficients 
that can be identified in each column. Emphasis is supplied on the 
word identified because the user generally must perform the identifying 
and numbering shown in Figures 2. 7-2. 8 and this may not be a trivial task. 
Another possibly troublesome feature of these solvers is that they usually 
do not assemble the element stiffness coefficients, again this nontrivial 
task is left to the user. 


(6 






9 








For entirely in-core column solvers the tasks of reduction and back- ^ 
substitution are straightforward and essentially the same as for band and 
frontal solvers. However, when the equations are too large to fit into 
core they must be blocked as illustrated in Figure 2,8. The reduction and 
backsubstitution phases here are considerably more complicated and unfor- 
tunately can require substantial I/O. In the simple illustrated example, 
the reduction of block 2 requires that only block 1 be available. Similarly, 

the reduction of block 3 requires that only block 2 be available. But 

block 4 needs blocks 1, 2 and 3! The minimum working storage for our 
example would be 2 blocks (about 40 words), and if this was the maximum 
available core, it would be necessary to write and read blocks 1 and 2 
twice. Thus, more than 1 word of I/O is required for each coefficient 
(e.g., as required in frontal schemes). 

Another obvious I/O inefficiency with column solvers is that the 

assembly and reduction phases are separated, thus requiring at least twice 

the I/O of frontal solvers. However, there are significant advantages to 
column solvers. For example, in the equations illustrated in Figures 2.6- 
2.7 only 58 words are necessary to store the active coefficients. Since 
the last column is full, frontal and band solvers would need to store the 
entire upper triangular region of 136 words. The column storage scheme 
would also require fewer computations for the same reason. This advantage 
is not effective in the regularly number cartesian grids generally used in 
2D or 3D continuum problems and for this reason column solvers may be more 
effective for special applications with sparse matrices or very general 
purpose computer codes like SAP [7] and ADINA [9]. 

To illustrate the comparison between column and frontal solvers 
consider the grid of Q8 elements shown in Figure 2.9a & b. Both grids 
have 4 elements across the minimum connectivity direction. The column 
numberings in Figure 2.9a are like a band solver would number, and Figure 
2.9b gives numberings similar to a frontal scheme. Figure 2.10 shows the 
skyline for the band numbering. In comparing storage in Figure 2.10 to 
the frontal scheme which would have a virtually constant front f»14, it is 
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seen that the number of terms in Ihe skyline is virtually the same as the 
number within the superimposed front. In Figure 2.11 a similar skyline is 
shown for the frontal numbering scheme of Figure 2.9b. Here there are about 
10% fewer terms in the skyline than the front. This potentially could lead 
to a 20% reduction in computations, but recall that at least twice as much 
I/O must be done by the column solver. 
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Figure 2.11 Skyline for Frontal Numbering 



CHAPTER 3 CHARACTERISTICS OF SPECIFIC SOLVERS EVALUATED 


3.0 General Remarks 


Due to the scope and funding of the present study, it was not pos- 
sible to obtain optimally coded solvers of all types. Since virtually all 
computer systems are unique and require specialized coding to operate at 
peak efficiency, this means that those imported solvers will be at a 
disadvantage to those that were coded on site. While it was not possible 
to extensively recode the imported routines, we did attempt to use these 
solvers fairly in our study. This does not mean that our final evaluations 
will be fair or impartial. In fortran programming there is always the 
not-invented-here syndrome to be contended with. 

Also, there was no attempt to make an exhaustive search to locate 
new and innovative routines. Rather, the author chose to use routines 
that were gene» .ily available and known to him a priori. Clearly, many of 
the better solvers may be not generally available because they are pro- 
prietary or because the developer does not code an easily transferrable 
modular routine [14]. 

The following is a brief description of the actual band, frontal and 
column solvers that were used in the present study. 

3.1 Bandsolvers 

3.1.1 BANSOL 

BANSOL is a general purpose out-of-core equation solver that stores 
and operates on the coefficients in bandwidth blocks as illustrated in 
Figure 3.1. BANSOL is a very general standard utility routine that was 
originally coded by E. L. Wilson and subsequently modified by the author 
when used in the early finite element codes PALOS [15] and NAOS [16]. It 
is an out-of-core solver in that the reduction takes place block by block 
with at most two blocks in core at a time. When the reduction of a block is 
complete, it is written onto low speed storage and the next block is formed 
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and reduced. The minimum storage required for Bansol is 2b^. If more 
sUrage is available the number of DOF that can be stored per block is 
increased according to the formula n^^ = NDIM/b - b where NDIM is the avail- 
able core storage and b is the bandwidth. 

For the grid shown in Figure 3.2, after element N is assembled, the 
equations up to point n^^ are complete. The coefficients that have been 
assembled up to element N are stored as shown in Figure 3.1 in two- 
dimensional array form using singlely subscripted arithmetic. Forward 
elimination for equations 1 to n^^ is done and then the revised coefficients 
from equations 1 to n^^ are written onto low speed disk storage and the 
remaining triangular portion is shifted upward so that equation n.+l in 
block 1 becomes equation 1 in block 2. The process is then continued as 
before assembling the stiffness matrix from the N+1 element. Note that 
BANSOL combines assembly and elimination and consequently avoids the 

extra I/O charges for the additional write and read required when these 
phases are done separately. 

BANSOL uses an efficient high speed binary I/O routine and does not 
use standard fortran I/O. However, since the usual block size is 10^ words 
or larger, this is not an important factor. 

BANSOL is generally useful for solving small to intermediate sized 
systems of equations associated with 2D elements without midside nodes 

(b£l50). When 2b exceeds the available core requirement, other solvers 
must be used. 

BANSOL was used in the study as a base line comparison for the inter- 
mediate sized grids and of course is not effective for large 2D or 3D problems 
with larger bandwidths. 

3.1.2 USOL 

USOL is an out-of-core band solver that has the capability of parti- 
tioning the band when the bandwidth blocks exceed the available core capacity 
USOL was developed by E. L. Wilson and was the original solver used in the 
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«Tba'Sol *" «*entiany the se™e manner , 

BANSOL. however, the process of assembly and reduction are separate No 

special coding was done on USOL to optimize its logic or I/O efficiency 

since It was anticipated at the outset of this research that bandsolvers that 

th,T, r '“P^titive. Also, it is generally known 

that USOL IS not an efficient solver for grids with midside nodes or for 

large bands. It was used in SAP because of its ability to use available core. 
3.2 Frontal Sol vers 


3.2.1 ZIPP 

ZIPP is the very general frontal solver developed by Irons [17] and 
use extensively in the TEXGAP computer programs [1]. In this report there 

^ed by Irons and listed in his paper [17]. The original version used Fortran 
/ nd was written to include many options for multiple load vectors and 
resolutions for nonlinear problems. The first modification converted all I/O 
and removed some burdensome logic related to the multiple load vectors and 
resolutions. This version is called HSZIPPl for High Speed binary I/O. 
ince ZIPP uses its own dimensioned buffer [17] there is greater control 
possible over when the internal buffer is released to the I/O buffer (or 

s^oIL'^ff ^ version released the dimen- 

of words t f ‘"“P ‘Pe tyPfcal number 

of words tr nsferred was 500-2000. A simple modification was made so that 

the internal dimensioned buffer was dumped only when it was full, thus 

5000 words are usually transferred. This second version is called HSZIPP2. 

As described in Section 2.2.2, ZIPP permits a very general node 
labeling scheme to be used because it has an extensive prefront routine 
that stores all DOF labels and then tests all other labels to determine 
the necessary information about the first, intermediate and last appearances 
a label. This is a computationally expensive procedure that increases 
as where N equals the number of degrees-of-freedom per element times 
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the number of elements (N=4n for 2D problems). This particular prefront 
feature is not common to all frontal solvers, for example, PUZZLE does not 
perform this extensive search. 

3.2.2 PUZZLE 

PUZZLE is a frontal solver developed by C.P, Johnson which can parti- 
tion the front when the front size f^/2 exceeds the available core [18]. 

Another important feature of this code is its extensive substructuring capa- 
bility. The frontal method intrinsicly generates substructure stiffnesses 
as the front progress, and the coding in PUZZLE fully exploits this useful 
feature. PUZZLE operates in a similar manner to ZIPP although the coding is 
quite different. The PUZZLE prefront is less general than ZIPP and is consi- 
derably more efficient. For typical continuum problems this loss in generality 
is minor. 

Like the ZIPP versions, PUZZLE handles all I/O using the high speed 
lOP routine and operates at near peak I/O efficiency on the Texas system. 

A very significant portion of this I/O efficiency occurs with large fronts 
that require partitioning of the front. Say for example that the front is 
almost constant and core memory is available only for about 1/2 of the needed 
space. Thus, 2 partitions are necessary for all frontal blocks. With USOL, 
COLSOL and most other band partitioning or column solvers this would at least 
twice as much I/O on the coefficient scratch file. However, PUZZLE is coded 
in such a way that needless repetition of I/O is minimized and frequently 
reduces the I/O burden to 1.3-1. 5 times that required If all coefficients 
fit in-core. For problems that require 3-10 frontal partitions, this savings 
could be 1-2 orders of magnitude. 

3.3 Column Solvers 
3.3.1 SKYSOL 

SKYSOL is an in-core skyline or column solver developed by C. A. 

Felippa [12]. Since SKYSOL is in-core, the maximum capacity is limited 
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to problems with 400 to 1000 DOF. In situations where there are only a 
few long columns, it is especially efficient in reducing storage as illu- 
strated in Figure 2.6. SKYSOL is useful for one-dimensional and some 
small two-dimensional problems. It is used in this research to provide 
''base-line'' comparisons, since it is anticipated that all out-of-core 
solvers will be less efficient. 

The method used for storing the coefficients is the compacted column 
storage scheme discussed in Section 2.3.2. All the elements above the 
main diagonal are stored in sequence by columns. All coefficients above 
the first nonzero coefficient in each column are omitted and all coeffi- 
cients between that point and the main diagonal are stored whether or not 
they are zeros. To use the compacted storage scheme, a locator array is 
needed to locate the position of the diagonal coefficients in the one- 
dimensional array. The locator array gives the necessary information to 
assemble the total stiffness and determines how many elements need to be 
revised during the elimination procedure of a given coefficient. 

3.3.2 GASP 

GASP is an active column or skyline solver coded by C. P. Johnson 
at the University of Texas at Austin primarily for the analysis of thin shell 
structures (which usually have 5 or 6 DOF per node). Since GASP was not 
intended for use in a general purpose code, the column heights are restricted 
to be such that they do not extend past the previous block. For example, 
in Figure 2.8 the column heights in block #3 extend only part way into 
block #2 but block #4 extneds back to block #1. In GASP, block #4 would be 
restricted to a height of 8. This restriction is important in limiting 
I/O since multiple passes through the stiffness coefficient scratch file 
are eliminated. 

GASP was coded using the high speed binary I/O routines available at 
Texas, thus further enhancing its I/O efficiency. As with all column solvers, 
GASP must separate the assembly and reduction phases and make multiple passes 
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through the element stiffness file. However, cumulative lower and upper 
limits on these passes are computed to reduce needless I/O. GASP requires 
the nodal connectivity array to be stored in-core and requires that the 
array contain DOF numbers for purposes of storage (gaps are permitted). 

3.3.3 COLSOL 

COLSOL is an active column solver coded by Wilson and Dovey at the Uni- 
versity of California to replace USOL in the SAP codes [8]. Like USOL, COLSOL 
can handle virtually unliii'ited bandwidths or column heights or conversely 
can fit into very small core when necessary. Thus, COLSOL must handle column 
heights that extend past the previous block, e.g., block #4 in Figure 2.8. 

All I/O in COLSOL is unformatted fortran I/O and because COLSOL was obtained 
late in this study it was not possible to improve it. For this reason, 
timings on COLSOL are likely to be slower than for more optimized coding. 
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CHAPTER 4 COMPARISONS BETWEEN THE VARIOUS SOLVERS 


4.1 Basis of Comparison 

The main thrust of this research is to compare the computational 
efficiency of the various methods of solving systems of equations by per- 
forming controlled numerical experiments on the solvers described in 
Chapter 3. The basis of comparison will be the "computer time" required 
by each solver to effect the solution of a set of model problems. 

It is extremely important in evaluating this particular research 
to at all times remember the nature of the reported computer timings. All 
results presented herein were obtained on the University of Texas 6600/6400 
operating system. This system permits two types of timing; a "central 
processor(CP) clock" and a charge time (TM). TM time equals CP time 
(sec) + 0.004 PRUs where 1 PRU =64 (60 bits) words read or written. All 
times are in seconds. Every possible effort was made to run the model 
problems on the 6600 because the 6400 gives "equivalent 6600" time. In 
practice, it is found that runs on the 6400 were generally 10-50% longer 
than on the 6600. In spite of our efforts to eliminate this variable, it 
is possible that some of the times reported are from 6400 runs. In any 
event, CDC 6600s are well known for their poor "clocks" and variations of 
+15% are possible depending on how many "rollouts" of the job are made. 
Generally, jobs are run during peak demand periods give longer times than 
those run at low demand periods. In spite of all these difficulties, it 
was still possible to adequately interpret and compare results and fre- 
quently duplicate runs gave variations as 1 >w as +1%, 

When comparing equation solvers it is important to establish a 
common basis of comparison, or the results of the experiments could be 
meaningless. To this end, all equation solvers were required to read a 
sequential file of element stiffness coefficients (including the right 
hand side), assemble the coefficients, reduce the coefficients (forward 
elimination) and backsubstitute to compute and print the solution. The 
actual element stiffness coefficients used were: 
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A-=* 1 .0, i-1 ,2,. . . .N 

A. j= 10'^°, i?«j=l,2,....N 

B. = 1.0, i=»l,2,....N 

This permitted a check on the solution since Xj=l to machine accuracy 
(14 figures) and avoided the special case of A^j=0 for i^j which could 
lead to erroneous comparisons since some solvers skip on zeros. The sequen- 
tial file of stiffness coefficients was generally written and read using 
standard FORTRAN unformatted I/O so that a minimum charge of 512 words 
(8 PRUs) was incurred for each read or write. Some of the runs (especially 
the ZIPP runs) were made using the high speed binary I/O routine lOP and 
this reduced these I/O changes by up to a factor of 10. However, since 
most of the problems required only one pass through the element tape this 
effect was negligible since most of the I/O takes place when the assembled 
coefficients are read and written. 

Equation solvers also need information on the connectivity of the ele- 
ments and in general this information was also passed to the solver by means 
of a sequential file. However, it was not possible to demand a like effort 
from each solver because of the particular coding in the solvers that were 
available. For example, ZIPP reads such a file and takes the information 
in the file as labels from which a frontal storage position is computed. 
PUZZLE performs a similar function but requires the labels to be integer 
numbers from 1 to N. Gaps are permitted, but storage is allocated for N 
numbers; a much more restrictive form than ZIPP which creates the numbering. 
USOL and COLSOL both require degree-of-freedora numbering a priori so that 
no bookkeeping need be done in the solver. Clearly, no meaningful compar- 
isons can be established for the connectivity-bookkeeping phase. However, 
this is not a troublesome point since in all solvers except ZIPP this time 
Is small. ZIPP does require a substantial amount of time in its prefront 
and for that reason it is reported separately. 
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Whenever possible, all identifiable times have been reported: these 
include: 

(1) prefront and/or bookkeeping 

(2) assembly 

(3) forward elimination 

(4) backsubstitution 

(5) total time In solver 

Because of their intrinsic nature, it is not possible to provide separate 
timings on each phase for each solver. For example, ZIPP and BANSOL com- 
bine assembly and forward elimination, but USOL and COLSOL separate them 
at the expense of an extra write and read of the entire system of equations. 

In spite of the aforementioned difficulties, it is still possible 
to make meaningful comparisons between the solvers by simply monitoring 
the total CP and TM time required for each solver. 


4.2 Comparison Problems 

The most Important Item In the comparison basis Is of course the set 
of sample problems. This research was directed at large out-of-core solvers 
and for that reason Intermediate and large- two-dimensional (2 DOF per node) 
and small three-dimensional (3 DOF per node) problems were chosen. It was 

not necessary to select Intermediate or large 3D problems as the results 
will show. 

At the present time the majority of finite element production prob- 
lems still seem to be solved with first order or linear elements, that Is, 
elements with only corner nodes, although use of the quadratic elements 
(elements with 1 midside node per edge) Is Increasing. For that reason, 
and because this research was partially aimed at the TEX6AP computer codes, 
the following element types were used: 


Q4-4 node quadrilaterals (8 DOF) 


Q8-8 node quadrilaterals (16 OOF) 
B8-8 node bricks (24 DOF) 

820-20 node bricks (60 DOF) 


cubir n rectangular grids In 20 and 

cable grids ,n 3D. These are highly specialised grids but they are probably 

representative of the majority of production problems. In any case it is 

doubtful that grid regularity affected the results appreclabl^. Unl’e 

fron s I of equations (e.g., extr^ely sparse 

Irsnl Identified and taken advantage of in 

olver. the results for Irregular grids will be similar. Of particular 
importance in the choice of model problems is the presence or absence of 
n.1d , e nodes. If there are only corner nodes in the grid (Q4 or B8 e e 
irents) then the front and band are equal. However, if midside nodes are 

less. Midside nodes are very troublesome to bandsolvers and run the band 

imately the square of the band. Hidside nodes cause no difficulty with 
frontal or column solvers both of which appear to treat these nodes effi- 
lently with frontal solvers perhaps having a small edge. Frontal solvers 

on*ra ' 7r - " nee 

on y a few DOF per element are eliminated. In general, band, frontal and 

Lerartr"^* r ‘=“"P“t»tlons for grids with no midside 

node so the solvers with the least bookkeeping will always perform best 

Bookkeeping is generally smallest in bandsolvers and highest in frontal 

4 1 andTp"? T' ‘hP Scids shown in Figures 

Test Prob em #1 was run on this grid with the numbering across the 10 
element direction so that the band(b) and front(f) widths were equal. 
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Band & Front for Prob #1 b»f»26 


Figure 4.1 Test Problems II & 2 for the Q4 
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b=f=26. Test Problem #2 was run on the same grid but numbering across 
the 24 element direction so that b=f=54. In Figure 4.2 is shown a grid 
of 5 X 12 = 60 Q8 elements and 11 x 25 - 60 = 215 nodes and 430 DOF Test 
Problem 43 was run on this grid with the numbering across the 5 element 
direction so that b=40 and f=32. Test Problem 44 numbers across the 12 
element direction so that b=82 and f=60. 

These four problems are not large but neither are they small in size 
and they were selected because the general resolution between the Q4 and 
Q8 grids will be comparable. Thus, for a given modeling the four problems 
give a range of bands and fronts while retaining approximately the same 
number of DOF. The missing 60 nodes in problems 3*4 could be handled by 
a frontal or column solver with little additional computational effort. 

4.3 Results 

The timings for problems 1-4 are given in Tables 4.1-4. Similar 
results are also given in Tables 4.5-6 for a 36 x 10 grid of Q4. and an 
18 X 5 grid of 08s. This had the effect of increasing the band and front; 
n=814, b=f=78 for Test Problem 45 and n=634, b=118, f=84 for Test Problem 46. 

The results reveal that as expected the in-core SKYSOL is the most 
efficient in both total CP * TM time for the solution. When considering 
TM time. BANSOL was always the next most efficient and usually was close on 
CP times. Note that BANSOL S USOL were faster for CP times than the frontal 
and column solver for problems 1. 2 S 5. the Q4 problems, but were slower 
for problems 3. 4 S 6, the Q8 problems. This is caused by the midside 
nodes running the bandwidth up and demonstrates the general superiority 
of frontal and column solvers fop grids with midside nodes. 

The CP times reported for USOL are generally comparable but the TM 
times are very poor. This is caused by high 1/0 changes from unnecessary 
I/O such as forming the equation separately before beginning the reduction 
process. This doubles the I/O changes because if sufficient core is avail- 
able to store the full band, front or active columns and no secondary parti- 
tioning is necessary, then each stiffness coefficient must be written and 
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SKYSOL 

BANSOL 

USOL 

ZIPP 

TM(sec) 

3.26 

6.94 

25.27 

61.08 

CP(sec) 

2.22 

3.05 

4.73 

7.26 

Working Array 
Storage 

13,441 

12,000 

12,000 

12,000 

I/O PRO'S 
TM(sec) 
0.004 sec/PRU 

3.54 

6.83 

29.34 

78.20 


HSZIPPl 

HSZIPP2 

PUZZLE 

GASP 

COLSOL 

40.98 

25.56 

26.91 



5.70 

5.08 

6.09 



12,000 

12,000 

12,000 



38.54 

23.72 

35,38 




Table 4,1 Test Problem #1 10x24 Grid of 04's 

n=550 b=f=26 



SKYSOL 

BANSOL 

USOL 

ZIPP 

TM(sec) 

5.56 

15.17 

37.77 

64.79 

CP (sec) 

4.53 

6.62 

4.36 

11.00 

Working Array 
Storage 

26,881 

12,000 

12,000 

12,000 

I/O PRO'S 
TM(sec) 
§.004 sec/PRU 

5.29 

10.48 

42.19 

78.20 


HSZIPPl 

HSZIPP2 

PUZZLE 

GASP 

COLSOL 

46.52 

29.60 

33.01 



10.25 

8.01 

10.50 



12,000 

12,000 

12,000 



39.50 

24.84 

28.16 




Table 4.2 Test Problem #2 24x10 Grid of Q4's 

n=550 b=f=54 






>' II fi III- . ii iiinn rillliiiitfiA— 'tia 



SKYSOL 

BANSOL 

USOL 

ZIPP 

TM(sec) 

3.28 

7.29 

13.81 

17,75 

CP(sec) 

2.22 

3.70 

3.95 

4.21 

Working Array 
Storage 

13,329 

12,000 

12,000 

12,000 

I/O PRO'S 
TM(sec) 
0.004 sec/PRU 

3.66 

6.46 

12.88 

20.57 


HSZIPPl 

HSZIPP2 

PUZZLE 

GASP 

COLSOL 

12.98 

9.46 

9.97 



3.26 

3.13 

3.08 



12,000 

12,000 

12,000 



11.77 

8.35 

13.53 




Table 4.3 Test Problem #3 5x12 Grid of Q8's 

n=430 b=40 f=32 


iialilMiHM 






i 



SKYSOL 

USOL 

ZIP? 

TH(sec) 

5.66 

34.00 

20.37 

CP(sec) 

4.63 

8.18 

6.54 

Working Array 
Storage 

25,929 

12,000 

12,000 

I/O PRO'S 
TM(sec) 
0.004 sec/PRU 

5,52 

28.83 

20.86 


HSZIPPl 

HSZIPP2 

PUZZLE 

GASP 

COLSOL 

17.37 

12.43 

13.01 



6.97 

5.04 

5.45 


- 

12,000 

12,000 

12,000 



12.40 

9.44 

12.03 




Table 4.4 Test Problem #4 12x5 Grid of Q8's 

n=430 b=82 f=60 



USOL 

ZIPP 

HSZIPPl 

HSZIPP2 

PUZZLE 

GASP 

COLSOL 

TM(sec) 

97,9 

103.4 

75.9 

54.4 

60.4 



CP(sec) 

14.8 

22.8 

19.8 

18.1 

24.4 



Working Array 
Storage 

12,000 

12,000 

12,000 

12,000 

12,000 



I/O PRO'S 
TM(sec) 
0.004 sec/PRU 

95.7 

116.6 

60.5 

40.7 

44.8 




Table 4.5 Test Problem #5 36x10 Grid of Q4's 

n=814 b=f=78 



USOL 

ZIPP 

HSZIPPl 

HSZIPP2 

PUZZLE 

GASP 

COLSOL 

TM(sec) 

83.97 

36.06 

28.92 

25.47 

24.99 



CP(sec) 

20.41 

14.02 

12.45 

11.92 

12.41 



Working Array 
Storage 

12,000 

12,000 

12,000 

12,000 

12,000 



I/O PRO'S 
TM(sec) 
0.004 sec/PRU 

67.53 

31.99 

19.02 

18.51 

18.40 




Table 4.6 

Test Problem #6 
n=634 

18x5 Grid 
b=118 

of Q8's 
f=84 




read only once. For this reason, no further studies on USOL were made since 
it is clearly inefficient. 

There are 3 columns of results given for the ZIPP solvers that reflect 
only improvements in the I/O handling and no computational changes. Note 
that reductions of up to 60% were made in the TM time and that this also 
resulted in CP reductions of up to 30%. Inefficient I/O is clearly to be 
avoided both to reduce wasted (I/O) charges and to increase CP efficiency. 

The results from HSZIPP2 and PUZZLE are seen to be similar as would be 

expected since both are frontal solvers and both use efficient high speed 
binary I/O. 

Another important variable reported in Tables 4.1-6 is the storage 
required for the solver with a nominal lower level cutoff of 12,000 words 
being taken. For Problem #1, SKYSOL needed only 13,441 words to store all 
coefficients entirely in core, but BANSOL or USOL would need n*b=14,300 and 
since only 12,000 were made available they needed at least 2 blocks of equa- 
tions to effect the solution. In Problem #2 almost twice the amount of core 
was needed by SKYSOL. BANSOL and USOL needed over 5,000 v/ords just for their 
2b blocks. ZIPP and PUZZLE actually needed only f^/2rl300 for their active 
storage. In Problem #4 SKYSOL needed 26,000 words; BANSOL and USOL needed 
2b =15500 just for the blocks, but ZIPP, PUZZLE and COLSOL could get by with 
only 3900 words. Problem US gets quite large and would have required 
2b =28,000 words for 2 full blocks and thus with only 12,000 available USOL 
had to effect a secondary partition - doubling its PRUs. Note that for 
Problem #6 USOL require over 3 times the TM time and 70% more CP time than 
ZIPP. ZIPP's storage requirements were still under 4,000 words, a factor of 
4 less than needed by USOL to stay in-core. 


4.4 Further Results for Frontal Solvers 
4.4.0 General Remarks 

The conclusion that can be drawn from the results given in Tables 
4.1-6 is that HSZIPP2 and PUZZLE consistently give faster CP and TM timings. 
For problems which do not require subdivision of the front there is little 
difference between HSZIPP2 and PUZZLE except for the more extensive prefront 
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used in the ZIPP codes. Because the frontal codes were found to be the 
more efficient, it was decided to study HSZIPP2 in more detail. Frontal 
solvers consist of 3 major parts; prefront, front and backsubsti tution. In 
the remainder of this chapter a study will be made of these parts. 

4.4.1 ZIPP Prefront 

The frontal technique in the ZIPP codes must perform extensive pre- 
front sorting on the connectivity to obtain equation numbers in the front. 
In so doing, it is necessary to sort through the connectivities element by 
element as shown by the symbolic Fortran coding in Table 4.7. For this 
reason, the CP times in the prefront should be uniquely determined by NIZZ. 
NIZZ = (No. DOF/element)*(No. of element). 

NLAST=0 

DO 300 NEL=1,NUMEL 
• • • • • 

DO 200 KL=1.KUREL 
NLAST=NLAST+1 
DO 100 N=NLAST+1,NIZZ 
IF (NIZlNLAST).NE.NIX(Nj) GO TO 100 
• • • • * 
loo CONTINUE 
200 CONTINUE 
• • • • • 

300 CONTINUE 

Table 4.7 Symbolic Fortran in ZIPP Prefront 

Since the sort is triangular the dependence should be proportional to NIZZ^. 
This is verified by the prefront timings shown in Figure 4.3 for the Q8 
elements. Although the logic given in Table 4.7 has some nominal dependence 
on KUREL and NELEM, the number of degrees-of- freedom per element and the 
total number of elements, respectively, the strong dependence is on NIZZ, 
the sum of all degrees-of-freedom in all elements. Thus, the CP times shown 
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2 4 6 8 10 

NIZZ(No. DOF/elm * No. elm) 


Figure 4.3 ZIPP Prefront Times for Ql 
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Figure 4.4 ZIPP Prefront Times for Va»"ious Elements 
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in Figure 4.4 for the Q4. Q12, B20 and L2 element types were not expected. 
Table 4.8 gives the relevant statistics for these element types (the ”L" 

type Is one-dimensional, “Q" types are two-dimensional and the ”B" type Is 
three-dimensional ). 


Element 

No. of Nodes 

No. of DOF 

Total No. of ; 
per element 



per node 

KUREL 

L2 

2 

1 

2 

Q4 

4 

2 

8 

Q8 

8 

2 

16 

Q12 

12 

2 

24 

B20 

20 

3 

60 


Table 4.8 Prefront Test ElemeMt Types 


As shown in Figure 4.4, for a fixed value of NIZ2, the L2 requires the 
most time, the Q8 and Q12 less time and the Q4 and B20 the least time These 
differences are substantial even for NIZZ as low as 4.000. There is no obvious 
explanation for these unusual groupings or the strong time differentials. 

The prefront times in ZIPP can be significant, but as discussed earlier 
this procedure is very general and may not be necessary for many problems. 

For example, for the B20 element with 3 degrees-of-freedom per node, it is 
necessary to do the sorting only on the 20 nodes rather than the 60 total 

degrees-of-freedom in the element. The node scheme will typically be 9 times 
more efficient. 

4.4.2 ZIPP Front 

Naturally, the majority of execution Is spent In the “front" portion 
of frontal solvers which consists of both the assembly and reduction phases. 
Front times have been measured for the Q4, Q8, 012 and 820 elements. 

In Figure 4.5 Is plotted CP time for the front portion vs. the total 
number of equatlons(n) for various front widths (f.60, 84,112 « 152) using 
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the Q8 element type. As would be expected, the time varies linearly with n » 
for fixed f. Figure 4,6 is a similar plot for the Q4 element type and reveals 
similar results. In Figure 4.7 is given a master plot of CP time per 1,000 
equations(n) vs. the front width(f) for fronts from 25 to 200 for the Q4, 

Q8, Q12 and 820 element types. Over this range in front widths, the curve is 
seen to mildly quadratic, but obviously for front widths between 50 and 200 the 
curve is nearly linear. This is further illustrated in Figure 4.8 where CP 
time vs. n is plotted. This indicates that ZIPP handles large fronts effi- 
ciently, perhaps improving in efficiency as the front grows (of course, the 
front must remain in core). It would be interesting to test further in the 
range of fronts from 200 to 500 (125,000 words of storage) to determine if 
this weak curvature continues. 

4.4.3 ZIPP Backsubstitution 

In frontal codes backsubstitution occurs element by element, and as 
in the prefront, backsubstitution times should be uniquely determined by 
NIZZ, with the variation nearly linearly. This is verified in Figure 4.9 
for all element types. Note that the vertical scale in Figure 4.9 reaches 
only 12 CP sec for NIZZ:14,000. This is substantially less than prefront 
times and can clearly be neglected. 

4.4.4 Further PUZZLE Results 

The results reported for PUZZLE show that it operates at an efficiency 
comparable to ZIPP for all of the test problems. Because PUZZLE is the only 
known frontal code to partition the front, further numerical studies were 
performed on a larger class test problem. A sequence of tests using the 

4x4x2 and 4x4x3 grids of 20 node brick elements (820, with 3 DOF per node) 

shown in Figure 4.10 were performed. The reason that 4x4x2 and 4x4x3 
grids were run was to determine the incremental time to solve a layer of 
4x4 820 elements. The front for these problems when numbering across the 
4x4 plane is 86 nodes or 258 DOF, thus the storage required for the entire 

front is 33411 words. The results for these runs are given in Table 4.9. 


SO 
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Figure 4.5 CP(sec) in ZIPP Front vs n for Q8 
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Figure 4.6 CP(sec) In ZIPP Front vs n for Q4 
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Figure 4.8 ZI 
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max/mi n 

max words 

TM Time (sec) 

no. partitions 

per partition 

4x4x3 

4x4x2 

A 

1/1 

33411 

130.8 

78.4 

52.4 

2/1 

33153 

187.1 

78.1 

109.0 

4/3 

14535 

289.1 

146.3 

142.8 

16/6 

8911 

314.3 

161.9 

152.4 


Table 4.9 Total TM Time in PUZZLE for B20 Problems 


The first line in Table 4.9 gives the TM times when sufficient core 
for the entire front is available and gives an incremental time per layer of 
4x4 B20's as 52.4 sec. The second line of the table gives the data for a 
similar run with virtually the same core, but slightly less than needed to 
fit the front entirely in-core. The column giving the maximum/minimum number 
of partitions reflects the differences in the number of DOF per element that 
are eliminated. As more DOF are eliminated, more partitions are needed. 

In general, the minimum number of partitions reflects the typical situation. 
The 4x4x2 grid with 2 partitions ran in essentially the same time as 1 par- 
tition because with only 2 layers the maximum front is not achieved. How- 
ever, the 4x4x3 required about 56 sec more TM time. Thus, the increase in 
A TM time was about 100%. The 3 partitions were run with about 15,000 words 
of core and the increase in A time over 1 partition was about 200%. However, 
the 6 partition run with less than 9000 words of core was almost the same as* 
the 3 partition run. Clearly, PUZZLE is very efficient at handling the extra 
I/O needed when partitioning the front. 

The typical prefront times for all runs (i.e., the rows in Table 4.9) 
was 10.2 sec for the 4x4x3, 6.5 sec for the 4x4x2 and 3.7 sec for the A. 

The typical backsubstitution time was 16.9 sec for the 4x4x3, 10.7 sec for 
the 4x4x2 and 6.2 sec for the A. 

A sequence of calculations were also run on Q8 grids with a variable 
number of elements across the front(n^). These results are given in 
Figure 4.11 where the ratio of PP/CP and TM time per element are plotted 
versus the total number of elements across the front. The front width is 
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easily computed from this by the formula 


f = 2n^+14 

Thus. Figure 4.11 gives results for fronts up to f=174. This requires only 
about 15,000 words of storage and consequently all runs were made without 
secondary partitioning of the front. 

PP stands for the standard CDC "peripheral processing" time. Note TM 
time equals CP time + a multiple of PP time. In PUZZLE the PP burden is quite 
low varying in an inverse fashion from 1.5 for n^=10 to about 0.25 for 
n^=80. TM time is seen to increase in a weak quadratic fashion with n.. 
varying from 0.20 sec/elm for n^=10 to 2.0 sec/elm for n^=80. This is^* 
very sim^ar to the behavior reported for ZIPP and again it would be very 

interesting to see if the weak quadratic behavior continues for larger 
fronts. 


ill Qualitative Evaluation of Column Solvers 

While the numerical experiments reported herein for the out-of-core 
column solver COLSOL were done as fairly as possible, these results are 
undoubtedly less than optimal since no effort was spent enhancing the I/O 
efficiency of COLSOL. For a more accurate comparison this must be accom- 
plished. However, there are some general qualitative conclusions that 
can be drawn. 

Column solvers generally require the assembly and reduction phases 
to be separated and this means that a minimum of 4 I/O transfers (2 reads 
and writes) must be performed. Thus, these column solvers start at a 
distinct disadvantage to frontal solvers like ZIPP or PUZZLE, Recall that 
for large fronts or bands the main consideration is I/O efficiency. Also, 
it is doubtful that central processor variations of more than ±10X will be 
observed between efficiently coded solvers of either type. Another dis- 
tinct disadvantage of column solvers is that for regular grids with large 
fronts or bands, the connectivity between column blocks will always be more 
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extensive than for frontal blocks because one frontal block will only con- % 
nect within itself when partitioned, but column connectivities may easily 
extend across one or more major column blocks. 

For these reasons it is doubtful that even the most I/O efficient 
column solver will ever achieve the low I/O changes reported herein for 
PUZZLE. This is not to say that column solvers are not useful for many 
applications where the band is highly variable, but the preliminary eval- 
uation is that for regular 2D and 3D grids (especially with midside nodes), 
frontal solvers will give lower charges. Notwithstanding, further study of 
column solvers is clearly indicated. 
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CHAPTER 5 CONCLUDING RFMAPK<; 


The results of this research are probably of questionable value 
because of the excessive I/O charges on the Texas computer and also because 
better control over the solvers is undoubtedly necessary. In the present 
study little attempt was made to recode any of the solvers (other than 
ZIPP) and the particular patches made to run solvers coded by outside 
sources are probably not optimal for purposes of careful and accurate com- 
parison. For this reason, the reader is cautioned to make similar studies 
on the particular compute system being used. However, it is the author's 
opinion that a great deal of useful qualitative information was learned in 
the present study (certainly by the author). 

Of relevant interest to those who solve large scale 2D and 3D 
continuum problems on a routine basis are out-of-core frontal and column 
(skyline) solvers. Clearly in-core and bandsolvers are useful only 
for testing or research purposes and have no place in large scale general 
purposes computer programs. There are two subclasses of these solvers of 
Interest, namely, those solvers that partition the front or active column 
area and those that don't (many column solvers automatically do this). 
Partitioning is required only when the front becomes too large to fitinto 
available memory. For example, if 4 x 10^ words can be dimensioned in 
central memory then a maximum front of 300 can be handled without parti- 
tioning. This would be sufficient for almost all 2D problems and some small 
3D problems. However, if for example large memory is used on a CDC 7600, 
up to 2.5 X 10 can be dimensioned permitting fronts of 700 to be processed. 
This would be sufficient for all but the largest 30 problems. Thus, users * 
at installations with large systems would typically use ZIPP which does not 
partition, but on smaller computer systems (e.g., Univ. of Texas, mini or 
micro-computers) PUZZLE or COLSOL would be necessary for large problems. 

The choice between frontal and column solvers has probably not been 
sufficiently resolved by the present study although the author has some 
preliminary evaluations. Selection between the two methods should currently 
be made on the basis of availability, specific adaptation to an individual 
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system (expecially the I/O), the nature of the code in which the solver is 
to be used and the particular bias of the user. The frontal method probably 
offers more flexibility to the programmer doing code development and in 
particular, ZIPP is quite easy to communicate with. Both frontal and column 
solvers offer strong potential for substructing and of course substructing 
invariably lowers computational costs and storage requirements when it can 
be used. PUZZLE appears to be a very strong competitor in the substructing 
area. At present the author's choice is the frontal solver, expecially for 
grids with midside nodes. 

The bottom line in evaluation of out-of-core frontal and column 
solvers appears to be quite simple. If partitioning is not done by the 
solver, then the I/O charges are fixed. There is no variation between good 
solvers of either type because the write and read the file of element stiff- 
ness coefficients and the file of assembled coefficients exactly once. A 
solver of this class that does more I/O is unnecessarily inefficient. There- 
fore, I/O efficiency is not a factor beyond the above point and efficiency 
can be measured solely on a CP basis. For solvers that partition the front 
or active column area almost the reverse is true because multiple passes 
through the element file and the scratch file for the reduced coefficients 
is required. Since this will typically increase I/O charges by a factor of 
2-10, I/O will almost always dominate the charge algorithm. CP efficiency 
cannot be ignored but variations of 25% would be in the noise level compared 
to I/O variations. PUZZLE appears to be an extremely strong competitor in 
the area of I/O efficiency for partitioning solvers. Unlike other schemes, 
PUZZLE does not require a integer multiple of passes through the scratch file 
and runs with multiples as low as 1.3 have been observed. 

The reader will have undoubtedly detected a bias towards CDC equip- 
ment and the author readily admits to such. The 60 bit word effectively 
eliminates concern over roundoff error for most problems. IBM equipment 
(32 bits) of course requires full double precision for all computations and 
Univac equipment (36 bits) likewise requires double precision for almost 
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all large systems of equations.* Double precision arithmatic is obviously 
more expensive and for that reason will alter any conclusions about CP 
efficiency. 

While CDC's word size is commendable » its I/O handling is certainly 
not. There are no standard I/O packages available on CDC equipment and 
it is this author’s experience that each CDC installation has either no 
such package or four of them; none adequately documented and generally total 
unknown quantities to the system programmers there. This is an unacceptable 
situation and at present effectively prohibits the transfer of I/O efficient 
codes. There are several important factors to consider when improving I/O 
efficiency. Most I/O routines do some form of buffering that may or may not 
be useful. Buffering is hardly ever useful and the author recommends direct 
transfer from central memory to disk without buffering. This is especially 
important when performing syncronous I/O to random storage (syncronous I/O 
waits until the transfer is complete before continuing execution). Random 
storage is generally desirable because the equation solver scratch file 
must be read in reverse order and because the block size does not need to 
vary, thus making direct addressing straightforward. Asyncronous I/O 
(execution continues while I/O is being done) appears attractive due to the 
obvious CP efficiency. However, in practice to effectively utilize this 
feature requires a fairly sophisticate coding effort and a tradeoff with 
the more straightforward syncronous transfer of single large blocks directly 
(without buffering) to disk may not give a clear advantage to asyncronous 
schemes. This entire process is future complicated by the CP interruptions 
on multi-processing computers. 

In closing it appears to the author that the question of out-of-core 
solvers is not a closed one. For static problems the choice seems clearly 
to go with a frontal or column solver. For implicit dynamic codes this 
choice is not obvious but the ’’tilt" would seem to go to iterative or hybrid 
methods, especially SOR and conjugate gradient techniques. Thus, there 
appears to be a need for further research. It would seem appropriate to 


♦The smallest recormended word size for our tasks is 48 bits. 
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continue the study of direct codes but develop frontal and column solver ' 
codes from scratch as opposed to using existing subroutines as was done in 
the present study. The author's choice would be to develop solvers like 
PUZZLE that partition when necessary but at no loss efficiency when parti- 
tioning is not required. A parallel effort on iterative methods also seems 
desirable at this time for several reasons. The next generation of computers 
will permit regular solution of 3D nonlinear dynamic problems. They will 
also be more prone to I/O binding than present computers simply because the 
CP will be faster. This will tilt the balance towards number crunching and 
away from I/O thus making iterative methods more attractive than at present. 
This latter task is clearly more difficult and the benefits may be slow in 
being realized. 
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APPENDIX A . 1 

BANSOL BANDED EQUATION SOLVER 


WTF.^<ON/tkl¥ BUNK 
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PROGRAM 8ANTE8T( INPUT, OUTPUT, TAPE2, TAPES) 
DIMENSION IY(0) ,8(256) .F( 16) 

DATA lOnF/2/,MDlH/12»00/,LOIM/86a/ 

COMMON 

1 MRAnO, NNP, ND, MMaX, NUMMLK, B(8b«), A(12P00) 

READ 5,NI,NJ,NP 
5 F0RMAT(315) 

M8ANDs(3*nU5)*ID0F 

numel=ni*nj 

NNPa(MOlM/M0AND*M3ANO)/2 

NUMNP3(2*NI + 1)*(NJ+1) + (NI4.1 )*NJ 

IF(NNP,GT,NUMNP) NNPsnUMNP 

NP03lD0F*NP 

NOsIOOF*NNP 

IF(ND,LT,MBAN0) stop 

NMAXsND*MBAND 

NDlM3l.OlM4NMAX<fMBAND<»MBANO 
PRINT 92,ND,NP0,NUMNP,NNP,NUMEL,MBAN0 
92 FORM/T(* 92 NO NPD NUMNP NNP NUMEL MBAND* , 61 1 3) 
DO 500 N=l, NUMEL 
IF(N.EQ.l) GO TO 3 
DO 90 Ksl,NP 
XlsN-l 
I2=I1/NI 

IFn2*NI,E0,Il) GO TO 3t 

IK32 

IF(K,E0,a,0R,K,E0,8) IKsl 
GO TO 39 
31 IK=NI+q 

lF(K,EO.<l,OR.K,Ea,B) IKs2*(Nl + l) + l 

39 IY(K)sIY(K)+IK 

90 continue 

GO TO 91 

3 READ 20, (IV(I),IaJ ,NP) 

20 F0RMAT(8I5) 

91 KsM 

00 10 1=1, NPD 
DO 10 Jsl,MPD 
KsK + 1 

S(K)=l,0E-B 
IF(I,EQ,J) S(K)=1.0 
10 FCI)=1.0 
NPS=NPO*NPD 
NCOUNT=NP*MPS+NPO 
call I0P(2Hwa,2,IY,NC0UNT) 

IF(N,GT,2) CO TO 530 

500 CONTINUE 

CALL SECOND(Tl) 
call J03INF0(8,J1) 

CALL block (NUMEL, N0IM,NC0UNT,I00F,NP, NUMNP) 
call J0HINF0(8, J2) 

CALL 5EC0Nd(T2) 

CP3FLOAT(J2-J1)/1000, 

DT»T2-TI 
PRINT 501,OT,CP 

501 FORMAT(* TOTAL TM TIME IN BANSOL » *,F10,5/ 

• * TOTAL CP TIME IN BANSOL = *,F10,5) 

end 

SUBROUTINE BLOCK ( NUMEL » NO I M,NC0UNT,I00F,NP, NUMNP) 
dimension UR(1000) 
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CQHHQN 

1 MBANO, NKP, NO# NMAX, NUMBI.K, B(664), 
DIMENSION LH(8)#IY(6),S(2S6)#Fn6) 

* FORM STIFFNESS MATRIX IN BLOCKS 


NUMBLKaB 
DO 50 NSl»NOlM 
50 B(N)«0,0 
ITY»l 
NEL30 

60 NUMBLKaNUMBLK-f I 
NMaNNP*NUM0LK 
NLsNMxNNPf 1 

KSHIFTsI00F*(NL-1) 

IF(NM,gE,NUMNP) iTYsa 

c 

select element in block 

c 

call I0P(3HREW,2) 

DO 210 Ns1,numEL 
CALL I0P(2HR8,2,IY»NC0UNT) 
IF(N,LE.NEL) CO TO 210 
DO 60 l3l,NP 

IF <iy(I),lt,nl) go To 80 

IF (IV(I),LE«NM) CO TO 90 
60 CONTINUE 
CO TO 210 
C 

ADD STIFFNESS AND FORCE VECTOR 

c 

90 DO 190 1=1, NP 

190 LM(i)= 2 *(iYcn*n 

NEL«N 

KKs0 

LLa0 

DO 200 I31#NP 

IIsLMCD^kSHIFT 

IJ3MBAN0*(II»1) 

NN=1J-LM(I)+1 
DO 200 K3i,2 
KKsKK^I 
II 3II^1 

IF(N,GT,5) CO TO 181 
161 0(II)3B(II)+F(KK) 

IJalJtMBAND 
NN=NN^MBANO»i 
DO 200 J=1,NP 
JJaNNfLM(J) 

DO 200 L81*2 
LLbLL+1 
JJsJJf 1 

IF (JJ.LE.IJ) go TO 200 

191 A(JJ)BA(JJ)fS(LU) 

200 continue 

210 CONTINUE 

NTsIOOF*NUMNP 
call 0ANSOL(UR,ITY) 

390 IF (ITY.NE,2) CO TO 60 
call BAN30L(UR#3) 

PRINT 391#(UR(I)#Ial,NT) 


A(12000) 
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wETUpy* uR(i)*,/nnEia,in 

End 

SUBHOUTInE H^NSOUdJ, ITY) 

c 

c....» I,n aANOfn EOIIATIO- SOLVta, reducer and BACKSUaSTITUTER 


COMMON 

NIJMBLK, 0(864), 
f?? I0TA(3^). IOTB(30) 

PRINT 111,NO,NUM0LK,ITY 
111 F0RMAT(* 11 I MO NUMBLK ITV*/,3I10) 

IF (ITY,EQ,3) go to 300 


C 

c 


REDUCE Equations and load vector 


A( 12000) 


N3|«mBaNO 

M50 

110 n=n+mbano 

MaM+ 1 

IF (N.GT.NMAX) GO TO 200 
DsA(N) 

IF (D,EO,0,0) GO TO 110 

IsN 

IJsN 

MJsM 

031, /□ 

FsB(M)*0 

B(M)3F 

00 130 L3?,MBAN0 

1 = 1 + 1 

U = IJ + M8AnO 

MJsMJ+1 

C3A(I) 

IF (C,EQ,0,0) GO TO 130 

8(MJ)sB(mj)«c*F 

CsC*0 

NKsI 


NJsIJ 

00 120 K=L#M0AND 
A (NJ) aA (NJ)«C*A (NK) 
NJsNJ+1 
120 NKaNK+1 
a(I)3c 

130 CONTINUE 
GO TO 110 

200 IF (ITy,EQ,2) RETURN 


C***** 

C 


WRITE REDUCED BLOCK OF EQUATIONS 


C 

c***** 

c 


IOsIoP(2hw0,3,0,nD) 

I0sl0P(2HWfl,3 , a,NMAX) 
I0=I0P(2HWR,3) 

shift Equations for next block 


NS0 

I=N0 
NJs0 
IJSNHAX 
210 NaM^l 
1 = 1 + 1 
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0(N)B8(n 

00 220 U=!»m04NO 
NJxNJ^I 
IJalJft 
A(NJ)aAdJ) 

220 A(IJ)30,0 

IF (n.lt.mbano) Go To 2ia 

N3N4-J 

00 230 iBNfNO 
B(I)a0,0 

00 230 1.31,m9aND 
NJ=NJ+1 
230 A(i>ij)s0,0 
RETURN 

***** BACK8U98TITUTE FOR UNKNOWNS 

300 NK3N0*NUHBLKfl 
NHOROSaND+NHAX 
N88=NH0R0S/6a+l 
N8aNSB*(NUM9LK*n 
3!0 NsND*l 
XJsNMAX^i 
320 Nsn* 1 
NKaNK«l 
l3N 

XJaIJ«M0ANO 

NJrlJ 

FsB(N) 

00 330 K32,M8AN0 

lal^l 

Nj8NJ+i 

330 Fsf»A(NJ)*0(I) 

B(N)3F 

U(NK)sF 

IF (N,GT,l) GO TO 320 

numblksnumhlk-i 

IF (NUMm.K,LE,0) RETURN 

***** shift LAST UNKNOWNS AND READ NEXT BLOCK 

XsNO 

Na0 

00 350 L*l»MBAND 

X»I*1 

Nsn^I 

350 B(I)3B(N) 

NSaNS«NSR 

l0al0P(2HSP#3»NS) 

IOslOP(2HR0,3,0,NO) 

I03I0P(2HRB,3,A,NMAX) 

IF(IO,EO,o) go to 310 
PRINT 351 

351 FORMAT!* ERROR ENCOUNTERED*) 

8T0P 

ENO 


USOL 


appendix a. 2 

BANDED EQUATION 




SOLVER 
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c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


TAPE8al00«,TAPE5»i00» 

, # T APEaa 1 000, T APE 73 1 000 , TAPE83 1300, TKPLOTb 1000, TAPE5aI<'^PUT. 

, TAPEbaoUTPUT) ' 


*4 


33 


«?]!-* **- ** ** ** ** ** ** ** ** ** ** ** ** ** 33 

8AP2 A STATIC ANALYSIS PROGRAM FOR THREE«DlMENSIONAL STRUCTURES 

REVISED march 1972 


33 33 33 33 33 33 33 33 33 33 33 33 33 

COMMON AC12000) 

COMMON /KKW/ 8B(2080) 

dimension T(7),F(lb),LCQNEC(lb),S(l6,l6),LM(l6) 

DATA LL/l/#ID0F/2/ 

data F/lb3l,0/ 


33 33 


33 


3 


3 


PROGRAM CAPACITY CONTROLLED BY THE FOLLOWING TWO STATEMENTS ,,, 

read 5 ,NI,NJ,NP 
5 F0RMAT(3I5) 

MTOTS12000 

IFLsLOCF(A) 

NUMNPa(23Nlfl)3(NJtl)+(Nl3l)3NJ 

NUMEL3NI3NJ 

NP232ANP 

MBAN03lD0F3(33Nl35) 

NEQ3NUMNP3IDOF 


INPUT JOINT DATA**IO ARRAY ON TAPE fl 


Nisi 

N23N1363NUMNP 

N33N23NUMNP 

N4aN3-fNiJMNP 

N5aN43NUMNP 

N6sN5<3NUMNP 

call PP(4LRFLP,IFL+MT0T) 

form element STIFFNESSES-*STIFF, on tape 2 "STRESS MATRIX ON TAPEl 

rewind I 
REWIND 2 
REWIND 4 

INPUT NODAL LOADS AND JOINT MASSES — " WRITE ON TAPE 3 

NEQBs(MTOT-a3LL)/(MBAND3LL+l)/2 
NBLOCK3(nEO- 1)/NEQB +l 
IF (NEQB.GT.NEQ) NEQBsNEO 
N3sN2<fNEQB3LL 
N43N3+b3LL 

WRITE (6,201) NEQ,MBAND,NE0B,N3L0CK 
FORM BB(K) MATRIX, FORM ELEMENT STIFFNESS— "ON TAPE 2 
ND323NP 

LROal3N03(NO'f 1) 

DO 100 Ns1,nUHEL 
IF(N,EO,n GO TO 91 
DO 90 Ksl,NP 
Il3N"l 
I2all/Nl 

IF(I23NI,E0.I1) go TO 31 
lKs2 
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1F(K,EQ,4,0R,K,E(3,6) IKol 
GO TO 34 
31 iKsNUa 

IFCK.EQ.y.OR.K.EO.B) lKsa*fNl + l)4-l 
39 LCONEC(K)sl.CONEC(K) + IK 

90 CONTINUE 
GO TO 92 

91 READ 2U# (LCONECCI)»Iel»NP) 

20 F0RMAT(8I5) 

92 DO 94 KslfNP 
Kla2*LC0NEC(K)-l 
K2a2*LCONEC(K^ 

KR132*K-1 

KR2»2*K 

LM(KRt)3Kl 

LM(KR2)3<2 

BB(Kl)3BB(Kn^F(KRn 

HB(K2)38B(K2)4>F(KR2) 

94 CONTINUE 

00 95 Ial,NP2 
00 95 J3I,NP2 
S(I»J)3l,ME-0 

lF(I,EO..n S(I»J)al«^ 

95 CONTINUE 

»*RITE(2) LRO,NO» (L^ ( I ) » 1 = 1 , NO) , ( (S ( I » J) » Jal , NO) , Isl , NO) 

100 CONTINUE 
C 

C FORM total stiffness MATRIX —ON TAPE U 

C 

NE2Ba2*NEQB 
N2sNltNE0B*HBAND 
N3 sn 2+NEQ8*LL 
N4sN3t4*LL 
NN2aNUNE2BAMBAN0 
NN3 sNN 2^NE28*LL 
NN4sNN3’f 4*LL 
call SECOND(T(l)) 

CALL J08lNF0(8,Jn , . 

call AOl)STF(A(Nl)#A(NM2)»NtJHEL»NBLOCK#NE29#LL»MBANO»Nl,NjfNP) 

call J0HlMF0(fl,J2) 

CALL SEC0Ni)(T(2)) 

C 

C SOLVE For displacement unknowns 

c 

NSB3(MBAN0fLL)*NEQ3 
N8B83NEa8*LL*(2'KM0ANO»l )/NEOB) 

IF(NS0B.LT,NSB) NSSBsnSB 

N43N3^NS0B 

CALL SEC'>ND(T(3)) 

CALL J0BINF0(8,J3) ^ 

call US0L(A(N1) »A(N3)r^I^^) , NEOB » HBAND» LL» NBLOCK »NSB» 4, 3» 7# 2# 2) 
call J0BINF0(8»J4) 
call 8EC0N0(T(4)) 

CSTF«FLOAT(J2-Jl)/l0i’R. 

CU8OL3FLOAT(J4.J5)/IO00, 

T8TFaT(2)-T(n 

TUS0L3T(«)-T(3) 

CPTOL»C8TF+CUSOL 
THTOL«TSTF+TUSOL 
PRINT 30»CSTF,CUSOL»T3TF,TU8OL 
30 FORMAT!* 30 CSTF CUSOL T8TF TUSOL *»4F10,6) 

PRINT 3S,CPT0L#TMT0L 
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35 F0flM4T(* TOTAL CP TIME IN USOL = *,F10.6/ 

• * TOTAL T>» TIME IN USOL s *»F1«S,6) 

*01 F0RMAT(3aH2 TOTAL NUMBER OF EQUATIONS 

1 /3«H PANOWIOTH 

2 /34H NUMBER OF EQUATIONS IN A BLOCK 

3 /3aH NUMBER OF BLOCKS 


«»X5, 

(# 15 , 


subroutine moves (IC ALL »NEOB,M,N0LnCK#LL»NE2B,U) 

COMMON /KKW/ BB(2000) 

DIMENSION U(NE2B»LL) 

NN3lCALL*NEQB 
DO 1 laUNEQB 
DO 1 L31*LL 
NPaNN^I 
U(I,L)aBB(NP) 

CONTINUE 

1F(M,EQ,NbLOCK) RETURN 
DO 2 l3l,NEOB 
DO 2 L3l»LL 
NQsnP+I 

U(NEQBfI,L)sBB(NQ) 

RETURN 

END 

SUBROUTINE ADDSTF (A,B#NUMEL»N0L0CK»NE2B,LL»MBAN0,NI,NJ,NP) 

FORMS GLOBAL EQUILIBRIUM EQUATIONS IN BLOCKS 

dimension A(nE2B,M8AN0)»SS(1 ) 

common /EM/ LR0fN0»L*i(2592) 

dimension 0(nE2B,LL) 

equivalence (SS#N0) 

NE0B3NE2B/2 

K3NEOB+1 

XsNBLOCK 

MHsSQRT(X) 

M03M8/2+1 

NEBB 3 MBANE 2 B 

MHsi 


NSHIFTsm 
REWIND a 


FORM EQUATIONS IN BLOCKS ( 2 BLOCKS AT A TIME) 

ICALLS 0 

DO 1000 Msi,nBLOCK ,2 
OO 100 I5WNE2B 
DO 100 Jsi,mbANO 
100 A(I,J)30, 

call B 0 VEB(ICALL»NEQa,M,N 8 L 0 CK,LL»NE 28 ,B) 

IF (M.EO.NQLOCK) CO TO 200 
ICALL»ICALL+2 
200 CONTINUE 
C 

rewind 7 

rewind 2 

NAs7 

NUmE«NUM7 

ND 82 *nP 

IF(MM,nE.1) go TO 75 

NA32 

NUmEbNUMEL 

NUH 730 
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75 00 NsI^NUME 

REAO(NA) LRO, (8S(I),Is1,lR0) 

DO 6(9frJ Isi^NO 
LMNsI«lm(I) 

IIsLM(n-NSHIFT 

IF (II.LE,M,DR,H,GT,NE2B) GO TO 600 
00 500 J*1,ND 
JJ«LM( J) tLMN 
IF(JJ) 500,500»360 
390 KK»lf!^0*J 

900 *Cn#JJ)«A(II,JJ)fSS(KK + I) 

500 CONTINUE 
600 CONTINUE 

determine if STIFFNESS IS TO BE PLACED ON TAPE 7 


C 


C 


1F(MM,6T,I) go to 700 
00 650 Is UNO 
IIaUM(I)«N8HIFT 
IF(lI,GT.NE2a,AN0,II.LE.NE8B) 
650 CONTINUE 
CO TO 700 

663 hRITE( 7) LR0,(SS(I),l3l,LRD) 
NUM7«NUM7fl 


GO TO 660 


700 continue 

WRITE (9) (CA(I,J),IsI,nE08),J=1,MBAN0), ((0(I,L)»IsunEQ0),LsUlL) 


1000 


IF(H,EQ,NBL0CK) go to 1000 

HMaMMtl 

nSHIFTsnSHIFT*nE29 


RETURN 

1002 FORMAT (4F10,0) 

2000 FORMAT(10m2STRUCTURE 12X 25HELEMENT LOAD MULTIPLIERS / 

• 10H LOAD CASE 9X IHA 9X IHR 9X IHC 9X IHD/) 

2002 format (I6»7X,4|Fi0.3) 
end 

^SUBROUTINE U80L (9,B,MAXB,NEOB,MB,LL,NBLOCK,NSB,NORG,NBKS.NTU 

dimension A(NSB),0(NSB),MAXB(NEOB) NT2,NRST) 


NCsMBfLL 

NBRs(M0«l)/N£OBfl 

INCsnEQB*! 

NMBsnEQB*MB 

N2SNT2 

NlsNTl 

rewind NORG 
rewind N8KS 


C 


c 


REDUCE EQUATIONS BLOCK-BV.BLOCK 


105 


DO 900 NsI^NBLOCk 

IF CN.GT.t .ANO.NBR.EO.n CO TO 110 
IF (NOR.EO.l) GO TO 105 
rewind N1 

REWIND N2 
NIsnI 

1F(N,EQ.1) nIsnORG 
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READ (Nl) A 
l|0 DO 300 laUNEQH 
OsA(I) 

IF(D) llSf3{10,120 
l!5 MaNEQB*(N«in+I 
write (b»ll 6 ) M,D 

116 FORMAT (33H0SET OF EOlJATIONS MAY BE SINGULAR / 

^ • 26H diagonal term of equation 16, 6 M EQUALS lPE12,a) 

120 llal 

DU 12S Js2,nC 
llall+NEUH 
125 A(II)aACII)/D 
C 

DO 130 jsI,NM 8 ,NE 08 
IF (A(J),NE,0.) MAXB(I)aJ 
130 CONTINUE 
C 

JL3U1 

IP (JL.GT.NEQB) go TO 300 
II = I 

DO 200 JaJL,NEQB 
IIsIItNEOB 

IFdl.GT.NHB) CO TO 200 
CaA(II) 

IF (C.EQ.U.M) GO TO 20tl 
CaC*A(I) 

C 

KKaJ-Il 

HAXaMAXBCI) 

DO 150 JJ3II,HAX,NEQ0 
150 A(JJ+KK)3A(JJ+kK)-C*A(JJ) 

c 

KK=J ♦NMB 
JJsI>fNMB 
DO 175 L=l,LL 
A(KK)3 a(kk)«C*A(JJ) 

KKrKK+NEQB 
175 JJsJJ>,vE!36 
200 CONTINUE 
300 CONTINUE 

WRITE (NBK 8 ) A,MAXB 

C 

SUBSTITUTE INTO REMAINING EQUATIONS 

C ■ 

DO 600 NNal,NBR 
1F(N4NN,6 T,nHL0CK) go TO 600 
NlaNl 

IF(N,EQ,1> NiaNORC 
IF(NN,EQ.N9R) NiaNORG 
READ (NI) B 
lLal^NN*NE08*NEQB 
DO 700 laUNEQB 
IlalL 

DO 6R0 KalfNEQB 

IF (II,GT,NMB) Go to 690 

CaA(II) 

IF (C,EO,0,0) GO TO 690 
CbC*A(K) 

HAXeHAXB(K) 

C 

KKal-ll 
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600 


DO 600 JJ=II,MaX,^£qb 
B(jJfKK)=H(JJ+KK)*C*A( JJ) 


650 


XK3I4NM8 

JJ«Kf NHB 

DO 650 L=t»LL 

B(KK)=8(KK)-C*A(JJ) 

KKsKKfNEOB 

JJsJJ+NEOH 


690 

700 


IIslI-INC 

lL=ILtN£Q8 


700 


750 

B00 


IF(JJ8R,nE, i ; 60 TO 750 
00 700 I=t,NSB 
An)BB(i) 

CO TO 800 
write (N3) 0 
CONTIiNJUE 


900 


HsNi 

N1SN2 

N2sm 


BACKSUBST I TUTION • RESULTS 


ON TAPE NRST 


905 


IS=LL*NEQ8 
NE8sNEO0*(nbR+1 ) 
NUMaM8R*N£aa 

MAXsnE8*LL 

DO 905 1=1, max 
a(i)s 0 , 

REwInO NRsT 

c-»— 

Do 1000 Ns1,n0LqC^ 
backspace N8KS 
read (N8KS) A,MAX8 
backspace NBKS 
00 910 l= 1,LL 
KaL*N£a 

00 910 J31,NJM 
I=K-NEQ8 
0(K)s8(I) 

KaK«l 


910 


920 


ISNMB 

DO 920 L = WU 
Ka(L»l )*NE8 
DO 920 J31 ,nEQB 
Iai + 1 

KsKtl 

8(K)sa(I) 


00 955 l3l,fjEO0 

J3NEQB^I»I 

MAXsmaX8(J) 

IF (A(J),EQ,f) ) 

00 950 L»l,LL 
KKa.T + (L«l)*NEB 
JJaKK+l 
ILsJ+NEQB 
C3B(KK) 


GO TO 955 
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DO iisIl#max,neq>j 
C = C-A(II)<»S(JJ) 

9a« JJsJJ-M 
95 a H(KK)3C 

955 continue 

c 

Iso 

DO 960 L=l,LL 
Ks (L-n *NKH 
DO 960 Jsi,NtQ9 
Ksk + 1 
lal + l 

960 A(I)3B(K) 

C 

WHITE (NR3T) (A(p,I = l,LS) 
PRINT 970, (A(I),Isi,LSl 
970 FORMAT!* 970 A ( I) * / , ( 1 0£ 1 2 , 3 ) ) 

1003 continue 

c— 

RETURN 

END 
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A'P P E N D I X B 

fortran listings of 

FRONTAL SOLVERS 
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APPENDIX B. 1 
FRONTAL EQUATION 


SOLVER 
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PROGRAM HSMAIN ( INPUT , OUTPUT , T APE 1 3 , T APE 1 « , T APE 1 5. T APE I fc r TAPE 1 9 , 
, TAPE5=1WPUT»TAPE6=0UTPUT) 

COMMON /GORIER/ NUMEL 

COMMON /TAPES/ NIN,N0UT»NTAPEI ,NTAPE2»NTAPE3»NTAPER»NTAPE5# 

'cOMMON/ZIP/LVAeL(2e)!KUREL!LPREO,LZ»tOEST(28),DUM(5),MVABt(l60O 

COMMON EU(25100) 

COMMON /lOlNFO/ JPRINT 
dimension LCONECC28) 


*** NFUNC(I»J) DEFINED *** 

NFUNCCI# t 


1 FORMATClHl///,* — NEW PROBLEM — — — 

15 FORMAT( ///5X»*NUMEL=*»I5#5X»*NEWRH3s*,I5#5X#*NELPAZs*»I10 ) 
20 FORMAT ( I5,/»(10I5) ) 

51 FORMAT(* 51 CONECTIVITV*, (815)5 
30 FORMAT ( 15,/, '5F10,05 ) 

JPRINT = 6 
NTAPE3 =13 
NTAPEa sia 
NTAPE5 =15 
NTAPE6 =16 
NTAPE9 =19 
NINs 5 
NOUT a 6 


1000 CONTINUE 

IO=H8IO{NTAPE3, 10,LVABL»0»0) 

XOcHSIO(NTAPER, 1O»ELPA,0,0) 

PRINT I 
READ 10,NI,NP 

READ 10» NUMEL# NEWRHS ,NELPAZ 
PRINT 15, NUMEL, NEWRHS ,NELPAZ 
IF ( NUMEL .LE, 0 ) STOP 
READ 20,KUREL» (LCONEC(I) , I=lf NP) 

NP2=2*NP 

LZ=KUREL*(KUREL+3)/2 

00 100 N s 1, NUMEL 

IF(N,EQ.l) GO TO 91 

00 90 Ksl,NP 

I1=N-1 

I2»I1/NI 

IFnZANi.EO.Il) GO TO 31 
IK=2 

IF(K,EQ,9,OR,K,EO,8) IK=1 
GO TO 39 
31 IK=NI+a 

lF(K,EO,a,OR,K,EO,0) IK=2*(NI+1 )♦! 

39 LCONEC(K)sLCONEC(K)+IK 

90 CONTINUE 
1F(N,GT,12) GO TO 91 
PRINT 51,(LC0NEC(K),K=1,NP) 

91 DO 95 J=l,NP 
Jl»2*J-l 
J2=2*J 

tVABL(Jn=2*LC0NEC(J)-l 

LVA0L(J2)=2*LCONEC(J) 

100 IO»H3ICfNTAPEO,l»LVABL»29,0) •'piAIlK 

DO 101 I=1*LZ 89 


EL(I)sl.0E«8 

101 CONTINUE 

00 102 ISUNP2 
NF=NFUNC(I»I) 

ELCNF)al,0 

102 CONTINUE 
LZKsUZ^KUREL 

DO 103 ZsLZKrLZ 
EL(I)=1.0 

103 CONTINUE 

00 200 NslfNUMEL 

200 I0=HSI0(NTaPE3,1,EL»LZ,1) 

IOsHSIO(NTAPE3,10,LVABL»0,0) 

IQsHSIO(NTAPE4I,10,EL>0»0) 

CAUL SECnNO(TK) 

CALL joeiNFO(e,ji) 

CALL ZIPP(NELPAZ) 

CALL JOBiNFOCe, J2) 

CALL SECOND (TW) 

CP=FLOAT(J2-*Jl)/1000, 

OTKWsTW-TK 
PRINT 201,OTKW,CP 

201 FORMAT!* 201 TM IN ZIPP = *»F10,6/5X,*CP IN ZIPP « *F10,6) 

60 TO 1000 

END 

SUBROUTINE ZIPP(NELPAZ) 

C 

C...,« FRONTAL SOLUTION METHOD BY BRUCE IRONS 
C 

COMMON /FORIER/ NUMEL 

common /tapes/ NIN,N0UT,NTAPEl,NTAPE2»NTAPE3,NTAPEa,NTAPE5» 

, NTAPE6,NTAPE7»NTAPE0,NTAPE9 

COMMON/ZIP/LVA0L(28),KUREL,LPREQ,LZ»LDEST(26),DUM'5),MVABL(160) 
COMMON /lOINFO/ JPRINT 
common ELPA(l) 

DIMENSION NIX(l) 

equivalence (NIX,ELPA) 

c 

NFUNC(I,J)3l+(J*(J^l))/2 
CALL SECOND(T0) 

CALL JO8INFO(8,J0) 

NBI631000 
NELZ3«62 
NIZZ s0 
MAXPAS0 
NVABZ30 
LCUREOS0 
NIXEND5NELPAZ 
LVEN0s28 
MVENO=160 
DO 232 Isl,MVENO 
232 MVABL(I)s0 

In»HSlO(NTApE6r 10«ELPA«0,0) 

IO»HSlO(NTAPE«r 10»ELPA,0,0) 

IO=HSlO(NTAPE5pl0fELPA,0,0) 

C 

C PUT ALL element NICKNAMES IN LONG VECTOR, NIX 
C 

JWHERE*5 

DO 10 NELEM a i, NUMEL 
lOsHSIOCNT 4PEfl,2,LVA0L»29,l) 

IF(NIZZ*KI,REL+NELEM,GT,NIXEN0) 60 

qn 


TO 130 


DO 8 Ib1,kUREL 
NIZZSN1ZZ41 
8 NlX(NIZZ)e..l.VA0L(I) 

NIX(NIXEND*l-NEUEM)aNlZZ 
in CONTINUE 
Nlal 

DO 26 NELEH=1,NUMEL % 

LPREQsLCUREQ 

LCUREQsNvABZ 

NZaNIX(NXXENO^l-NELEM) 

KURELaNZ-NUl 

LZ«NFUNC(2*KUREL,KUREL) 

C 

C FIND NEW NICKNAMES AND USE EXISTING DESTINATIONS, 

c 

DO 22 NEWsNlfNZ 
NICsNIX(NEW) 

LOESaNIC 

IFCNIC.GT.B) GO TO 20 
DO 14 LDES=1,MVEND 
IF(MVABL(LOES),EQ,0) GO TO 16 
14 CONTINUE 
JWHERE=7 
GO TO 130 

16 MVABL(LOES)sNIC 

IFCLOES.GT.MAXPA) MAXPAaLDES 
KOUNTal 
C 

C RECORD FIRST, LAST AND INTERMEDIATE APPEARANCES, 

C 

C 

C***** THE NEXT FIVE(5) STATEMENTS SHOULD BE REPLACED WHEN 
C***** converting the COMPASS ROUTINE FIRLAS TO FORTRAN 
C 

NNOaMZZ-NEW*l 

NiCDaNiC 

LDESOaLOES 

CALL FIRLAS1:nnO,NICD,NIX(NEW),LDES0»LAST, MOUNT) 

LASTaLAST^NEW-l 

C 

KOUNTbkOUnT«NBIG 
NIX(LAST)aLDES * NBIG 
LDESaLDES^KOUNT 
NlX(NEW)aLDES 
20 LOESTtNEW-NUDaLOES 
22 CONTINUE 
NlaNZ+1 
C 

C RECONSTRUCT ELEMENT NICKNAMES AND COUPLE WITH DESTINATION VECTORS. 

C AND INITIAL ELEMENT STIFFNESS AND LOAD DATA. 

C 

DO 24 KLal»KUREL 

CALL COOEST(LDES,LDEST(KL)fNSTRES,NBIG) 

LVABL(KL)a.MVABL(LDES) 

IF(N8TRE8,NE,0,AND,NSTRES,NE,1) GO TO 24 
MVABL(LDES)30 
NVABZaNVABZ^l 
24 CONTINUE 
C 

C *** REWRITE all element INFORMATION ON TAPE 
C 

I0»H810(NTAPE5, l,LVABL»64,0) 

91 
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26 CONTINUE 

I0SHSI0(NTAPE5.9) 

CALL J08lNF0(fl,Jl) 

CALL SECONO(Tl) 

CP3FLOAT(JI»J0)/1000, 

OT=T1»T0 

WRITE(NOUT,870) DT,CP»NIZZ#NVAB2»MAXPA 
87P FORHAT(*0TIME IN PREFRONT s *F8,3// 

• * CP IN PREFRONT s *F8,3// 

, A NIZZ s *15// 

, * NVABZ 3 *15// 

, * HAXPA s *15) 

IO=H3IO(NTAPE5,10,ELPA,0,0) 
l05hSI0(NTAPE3, 10,ELPA,0,0) 

PRE^PROGRAM ended and element TAPE WRITTEN, 

ESTABLISH STORAGE REQUIREMENTS AND AREA BOUNDARIES 

NPARsnFUNC(0»MAXPA*1)*NELZ 

NPAZsNPARfMAXPA 

NBAXOsNPAZ*! 

NBUFFA3NELPAZ*NBAX0 

JWHEREsq 

IF(N0UFFA,LT,MAXPA+a) GO TO 130 
NRUN03NPAZ-MAXPA 
IBAsNBAXO 
DO 38 1s 1,NELPAZ 
ELPA(I)s0,0 
38 CONTINUE 

KURPAa0 


SEEK and assemble NEW ELEMENT, 


(10 


az 


4(1 

46 


IPOS30 

NWORDS30 

DO 92 NELEM=1,nUmEL 

I0=H8l0(NTAPE5,2,LVABL#64,n 

I08HSI0(NTAPE3,2»ELPA,LZ,1) 

L30 

DO 40 KL3l,KUREL 

CALL C00EST(L0ES,LDEST(KL) »NSTRES»NBIG) 
MVABL(L0ES)3LVA6L(KL) 

LVabl<KL)=LDES 

IF(LDES,Gt,KURPA) KURPAsLDES 

continue 

DO 64 LHSRHSs1,2 
LHS52I.LHSRHS 
1RHS3ULHS 
MNM*LH8*KUREL*1RHS 
DO 64 KL3l,MNH 
GO TO (42r44) fLMSRHS 

kgslvablckl) 

MGOaNFUNC(0,KG)*NELZ 
GO TO 46 

MGOsCKL*! )*HAXPA*NPAR 
HMN3LM8*KL*IRHS*KUREL 
CO 64 lL3i»MMN 
IOslVABl(IL) 

LsLaJ 

CE»ELPA(L) 

MG3MG04IG 

q? ■ 
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CONTINUE 

C 

C elihin.te variable in position loes, ano write equation for tape 
00 90 KLst»KUREL 

IF(NSTRES,NE,0,AND,NSTRES,NE,n GO TO 90 

NDEQNalBAfKURPA^a * ^ “ 

IFCNOEQN.UE.NELPAZ) go to 70 
EtPAClBA^DsNWORDS 
NWOROS2IBA-N0AXO+2 
NSsNkORDS/69 

1F(6«*NS,NE,NW0R0S) NSsNS^l 
IPOSsIPOS+NS 

3002 ^o^^^?{^^;S?ERiESJi^E'??PE;^SS^^^ 

J^’j;5I0(NTAPE6,l,ELPA(NBAX0),6a*NS,l) * 

NDEQNaI6A^KURPA-f4 
70 iBOlAOalBA^LOES 

ND1AG3NFUNCC0,LOES+1) ♦NELZ 
PIVOTsELPACNOIAG) 

ELPA(NOIAC)a 0,0 

JWHERE313 

NlCsMvABiaOES) 

1F(P1VOT,EO,0, J GO TO 130 

HGZsNEUZ 

JGZsKURPA 

IBOsIBA 

DO e« LHSRH3ai,2 
IF(tHSRHS,E0,2) JGZ=1 
DO 8fl JGai,jGZ 
IBAalBAtl 

CO TO (72,76) »LMSRH3 
72 MG03HGZ 
MGZSMGOtJG 

IF(L0ES,GT, JG) GO TO 7A 

MGrMGO+LOES 

GO TO 78 

7« HGaNFUNC(JG,LOES)+NELZ 

GO TO 70 

76 PGOa ( JG«1 } *MAXPA4NPAR 
MG3MG0+LDES 
MGZaHGOfKURPA 
70 NOELTalBO-MGO 
CONSTaELPA(MG) 

ELPA (IBA)aCONST 
IF(CONST,EO,0) CO TO BA 
CONSTaCONST/PrVOT 

ELPA(MG)30,0 

IF (UHSRHS.NE.t) GO TO 60 

simultaneously, create a simple rouno off criterion, 

60 NNLaHCOfl 

**** the call inner should be replaced when 
**** converting from COMPASS TO FORTRAN 

CALL I^NER(MG0,HGZ, CONST, ELPA(NNL),N0ELT) 

‘93 


ea CONTINUE 

ELPA(IBOUG)aPIVOT 

XBA3NDEQN 

ELPA(XBA)=KURPA 

ELPAdBA^DsLOeS 

ELPA(IBA-2)sMVABLaDES) 

66 MVABU(LDES)s0 

IF(MVA8L(KURPA) ,NE,0) GO TO 90 
KURPAsKURPA-1 
IF(KURPA,NE,0) GO TO 88 
90 CONTINUE 
92 CONTINUE 

I0aHSI0(NTAPE6,9) 

I0sHSI0(NTAPE6, 10) 

CALL J0BINF0(8, J2) 

CALL SECONDCT2) 

CPsFLOAT(J2»Jl)/l000, 

OTST2-T1 

HRITE(NOUT#87l) OT , CP » IPOS » NPAZ 
871 FORMAT(»0TIME IN FRONT = *F8,3// 

, * CP IN FRONT = aF8,3// 

, * NUMBER OF SECTORS a *15// 

, * NPAZ = *15) 

C 

C BACK-SUBSTITUTE INTO EQUATIONS# TAKEN IN REVERSE ORDER 
C 

NBZalBA 
NEQsNVABZ 
LPREOaLCUREQ 
NELEMsNUMEL 
lOsHSIO(NTAPE9, 10) 

100 IF(1BA,NE,NBAX0) go TO 102 
NSsNW0R0S/6il 

1F(6<1*NS,NE,NW0R0S) NSaNS+l 

IPOSalPOS-NS 

NBZaNW0RDS4-NBAX0-2 

IO=HSIO(NTAPE6,6#ELPA#IPOS,0) 

IOaHSlO(NTAPE6,2#ELPA(NBAXO),6U*NS#l) 

NWORDSaELPA(NBZ*l) 

XBAaNBZ 

102 KURPAsELPA(IBA) 

L0E53ELPA(IBA-1) 

NICsELPA(IBA-2) 

XBARsIBA»a 

XBAalBAR-KURPA 

XBDIAGsXBATLOES 

PIVOTsELPA(IBOIAG) 

ELPA(IBOIAC)80,0 

MGOaNRUNO 

HGZsHGO^KURPA 

CONSTsELPA(XBAR*1) 

NOELTsXBA-MgO 

NNL»MG0*1 

C 

C****« THE call 0SUB SHOULD BE REPLACED WHEN 
C***** CONVERTING FROM COMPASS TO FORTRAN 
C 

CALL BSUB ( NNL » MGZ# CONST »ELPACNNL)#N0ELT) 

C 

C place ANSWERS AND PREPARE FOR NEW ITERATIVE LOOP, 

C 

ELPA(MGO*LOES)aCON8T/PXVOT 


ono oooooo 


ELPA(IBDIAG)sPIVOT 

NEQsnE0«1 

108 IF(NEQ,NE,LPRE0) go to 100 

IO=HSIO(NTAPE5,6,LVARL,NELEM-1,0) 

I0sHSI0(NTAPE5,2,LVABL,6a,t) 

00 112 KL=1»KUPEL 

CALL CODE8T(LDES,LDEST(KL),NSTRES,N0IG) 

MPUNsNRIJNO+LOES 

ELPA(KL)sELPA(NRUN) 

112 CONTINUE 

ELPA(6«)skurEL 

I0«HSI0(NTAPE9, 1,ELPA(1),60,1) 

NELEM3NELEM-1 
IF(NELEM,NE,0) go TO 108 
I0=HSI0(NTAPE9,9) 

CALL J0BINF0(8, J3) 

CALL 5EC0N0(T3) 

CP=FLOAT(J3-J2)/1000, 

CP3aFLOAT(J3»J0)/1000- 

DT3T3-T2 

T3ST3-T0 

WRITEtN0UT,872) DT,CP,T3,CP3 
872 FORMAT(*0TIME IN BSUB a *F8,3// 

• * CP IN BSUB s *F8,3// 

• * TOTAL TIME IN ZIPP 3 *F8,3// 

• * TOTAL CP IN ZIPP 3 *F8,3) 

RETURN 

130 CONTINUE 

WRITE (N0UT,83«) JWHERE,NIC,PIVOT,LOE8,LZ,NELZ,NELEM, 

1 nbuffa,ni2z,nelpaz,maxpa,npaz 

CALL error (J WHERE » 0, 0, 10HZIPPJWHERE) 

83a FORMAT(9H0 JhHERE s,I 3,5X,5HNIC s,i8,5x, 

I =.IS,a*,6HMEU =,IS/ 

3 6HNIZZ s,l5,qx,0HNELPA2 =,I5,ax, 
a 9H MAXPA s ,ia,10H NPAZ s , 17 ) 

END 

SUBROUTINE COOE8T(LDES#LOE8T,NSTRE8#N0IG) 

NSsLOEST/NBIG 

LDES8L0EST*NS*NBI6 

NSTRES3NS-1 

return 

END 

function HSIO(NTAPE,NOP,A,NWOROS,NWAIT) 

••ft COC 6600/6900 HIGH SPEED BINARY I/O FORTRAN INTERFACE ROUTTwp 
KWAIT option not operational on COC, BUT IT IS USEFUL ON 
.... UNIVAC 1108 SERIES NITH THEIR NTRAN ROUTINE 

NWAITaO FOR NO WAIT I NWAITsl FOR WAIT To'cOHPLETE I/O 


C0MM0N/TAPE8/NIN,N0UT 
INTEGER H8I0 
IF(N0P-2) 1,2,3 

WRITE NWOROS STARTING AT ADDRESS A 


1 HSI0s10P(2HWB,NTAPE, A, NWOROS) 
GO TO 90 


READ NWOROS STARTING FROM ADDRESS A 

95 


c 


2 H3IO»IOP(2HR8,NTAPE,A,NWOROS) 
lF(HSIO,EO,e) GO TO 90 

C 

C,.... ERROR 0IAG0N8TIC FOR INCOMPLETE I/O TRANSFER 
C 

WRITE(NOUTp3000) NTAPE,N0P,NW0R08»NWAIT,HSI0 
3000 FORMAT(*0H8IO READ ERROR / NTAPE NOP NWORDS NWAIT H8I0*// 

• *615) 

CALL ERROR(2,0,0,4HHSIO) 

3 IF(N0P«9) 6,9,10 
C 

C..,,, SET POSITION TO SECTOR NUMBER NNOROS • USED FOR BACKSPACE 

w 

6 HSI0«I0P(2HSP,NTAPE,NW0RDS) 

GO TO 90 
C 

C.,.,, WRITE EN0*»0F^FILE MARK 

c 

9 HSlOaIOP(2HWF, NTAPE) 

GO TO 90 
C 

C.,,,, REWIND 

C 

10 H8IO«IOP(3HREW, NTAPE) 

90 RETURN 
END 


INNER 

* 

* 

* 


IDENT Inner 
entry inner 

BSSZ 1 


(MGO, MCZ, CONST, ELPA,NDELT,) 


COMPASS VERSION OF THE FOLLOWING LOOP 


* DO 62 1»NNL,MGZ 

* ELPA(I)sElPA(I) - CON8TaELPA(I+NOF.LT) 

* 82 CONTINUE 

* 

* INITIALIZE 


SB7 

1 

Sa3 

B3 

SAl 

61 

SBl 

XUB7 

SA2 

B2 

SB3 

X2 

SA4 

64 

Sa5 

B5 

SB5 

X5 

SAS 

A4*B5 

RX7 

X3*X5 

BX0 

X4 

GE 

B1,B3 


* inner loop 

* 


LOOP 


FX6 X0-K7 

NX6 X6 

8A6 AO 

8A4 A4fB7 

8A5 AS>B7 

RX7 X3*X5 

BX0 X4 


STORE 1 INTO B7 
CONST TC X3 
MGO TO XI 
MGO TO B1 
MGZ TO X2 
MGZ TO B3 
ELPAIMGO^n TO X« 

NDELT TO X5 
NDELT TO B5 
ELPAtl+NDELT) TO X5 
FIRST MULT, C0NST*ELPA(I+NELT) 
COPY X« INTO X0 
SKIP LOOP IF HGO^laHGZ 


ELPAfCON8T*ELPA(ItNDELT) 
normalize RESULT IN X6 
STORE RESULT IN ELPA(I) 
NEXT ELPA(I) TO X4 
NF.XT ELPA(IfNOELT) TO X5 
CONST* ELPA(IfNDELT) 
ELPACI) TO X4 


SBl BUB7 

LT BWB3fL00P 

* 

* LAST SUBSTRACTION 

* 

exit FX6 XB-X7 

NX6 X6 

$A6 A4 

EO B0,B0, INNER 

END 

XOENT BSUB 

entry BSUB 

BSUB BSSZ 1 


ADD I TO INCR COUNTER 
TEST FOR END OF LOOP 


normalize X6 

STORE LAST RESULT IN ELPA(I) 
RETURN TO PROG, 

(NNL,MGZ, CONST, ELPA(NNL),N0ELT) 


compass VERSION OF THE FOLLOWING LOOP 
DO 10« IaNNL#MGZ 

CONSTsCONST - ELPA(I)*ELPACI+NDELT) 
10 a CONTINUE 


* 

* 


* 

LOOP 


* 

EXIT 


INITIALIZE 

S87 

1 

SA3 

B3 

SAl 

B1 

SBl 

XI 

SA2 

B2 

SB2 

X2 

SAO 

B4 

SA5 

65 

SBS 

X5 

SA5 

A44-B5 

RX7 

X4*X5 

BX0 

X3 

GE 

Bl, 02, EXIT 

INNER LOOP 

FX6 

X0»X7 

NX6 

X6 

s;4 

A4^07 

SA5 

AS437 

RX7 

xa*X5 

BX0 

X6 

SBl 

flUB7 

LT 

Bl,B2,LOOP 

LAST SUBSTRACTION 

FX6 

X0^X7 

NX6 

X6 

8A6 

A3 

EO 

60, B0, BSUB 

END 

IDENT 

firlas 

ENTRY 

fxrlas 


STORE I INTO B7 
CONST TO X3 
NNL to XI 

nnl to B1 

MCZ TO X2 
MCZ TO B2 
ELPA(I) TO Xfl 
NOELT TO X5 
NDELT TO B5 
ELPA(UNDELT) to X5 

first multiply 

COY CONST TO X0 
SKIP LOOP IF NNL sMGZ 


C0NST.ELPA(I)*ELPA(I4NDELT) 
normalize result in X6 

NEW ELPA TO xa 
NEW ELPA(I4N0ELT) TO X5 
ELPA(I)*ELPA(ItNDEL) 

CONST TO XO X0 
ADDS ONE TO COUNTER 
TEST FOR END OF LOOP 

LAST SUBXTRACTION 
NORMALIZE RESULT IN X6 
STORE RESULT IN CONST 
RETURN 

(NNO,NICO, NIX (NEW), L0E8, LAST, KOUNT 


COMPASS VERSION OF THE FOL, LOOP 

KOUNTsl '-.'IGtNAL PAGE IS 

DO la LA8 bNEW,NIZZ of POr- 'y^jrTY 

IF(NIXClA8),NE,NIC) CO TO 16 

NIX(LA8)»L0E8 

LASTbLAS 
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J0 


KOUNTcKOUNT ♦ I 
continue 


* 

« 


* 

initialization 


firlas 

BSSZ 

1 



SAl 

Bi 

NND TO XI 


set 

XI 

NNO TO Bl 


837 

1 

1 TO 87 


SA2 

B2 

NIC TO X2 


8A3 

B3 

NIX(l) TO X3 


SAa 

B4 

LOES TO xa 


BX6 

Xa 

COPY LDES in X6 


SB2 

1 

SET B2al 

if 

8B4 

B2 

KOUNT»l 

LOOP 

1X7 

X3>X2 

NIX(LAS)-NIC TO X7 


NZ 

X7, INCR 

IF NIX(LAS),NE,NIC JUMP 


SA6 

A3 



SB2 

B7 

LASTsLAS 


SBa 

BAtl 


INCR 

S07 

B7»l 

INCREMENT LOOP COUNT BY 1 


SA3 

A3^1 

NEXT NIX (LAS) 

n 

LE 

B7, 01, LOOP 

if(las,le,nnd) loop 


SX6 

B2 

B5 TO X6 


SA6 

B5 

LAST TO B5 


SX7 

B<1 

KOUNT TO B6 

if 

SA7 

B6 



EQ 

B0,B0, FIRLAS 

RETURN 

A 1 1 M 

END 




SUBROUTINE ERROR (NERR* WORD) 

COMMON /TAPES/ NIN#N0UT 
IF(I,NE,0) go to 5 
WRITE(NOUT»10) NERR,wOR0 
STOP I 

5 WRITECNOUT# 10) NERR»WORD,I» J 
10 FORMAT(* ERROR NUMBER « *, 16,2X, A10,216) 
STOP I 


APPENDIX B.2 

PUZZLE FRONTAL EQUATION SOLVER 
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PROGRAM PUZZLE ( INPUT » OUTPUT , TAPE I # TAPES, T APES , TAPEa , TAPES , TAPE6 ) 
COMMON NDOF(2SO0) ,IQ(20),ME,PE(60),ST(60,60) 

COMMON / SUBT / N0EL,N0PT,ISAV,NSTR,MSUB,NUEL,NSUB,ISUB(63) 
COMMON / 0000 i MAXF,MAXB,MAXW,1PRNT 
COMMON / TRAC / NWS,NWR,NHS,NUP 

CALL second ( Ti ) 

CALL J08lNF0(8,jn 

CALL SETDATA 

CALL J0BINF0(6,J2) 

CALL SECOND ( T2 ) 

CALL FRONTIO 

CALL J0BINF0(8,J3) 

CALL SECOND t T3 ) 

CALL FRONTST 

CALL J0BXNF0(8,J4) 

CALL SECOND ( Ta ) 

CALL BACKSUB 

CALL J0BINF0(8,J5) 

CALL SECOND ( T5 ) 

TO = T2 - Tl 
CPsFLOAT(J2-Jl)/l000, 

PRINT 100, TO, CP 
TO a T3 - T2 
CP=FLOATtJ3-J2)/l000, 

PRINT 200, TO, CP 
TO a T<l - T3 
CPaFLOAT(J«-J3)/l000, 
print 300, to, CP 
TO a TS - Tfl 


CPaFLOAT(JS-ja)/l000, 
PRINT ^00,TQ,CP 
TO a T5 - Tl 
CPaFLOAT(J5-Ji)/1000, 
PRINT 500, TO, CP 


PRINT 600,NWS,NWR 
PRINT 700,NBS,NUP 


100 FORMAT 
200*FORMAT 
300*FORMAT 
<I00*FORMAT 
500*FORHAT 



*1* 5X 

*TIME 

FOR 


6X 

* CP 

FOR 


6X 

*TIME 

FOR 


6X 

* CP 

FOR 


6X 

*TIME 

FOR 


6X 

* CP 

FOR 


6X 

*TIME 

FOR 


6X 

* CP 

FOR 


6X 

*TIME 

FOR 


SETDATA* F13,3 / 
SETDATA* F13.3) 
FRONTIO* F13,3 / 
FRONTIO* F13*3 ) 
FRONTST* F13,3 / 
FR0NT8TA FI3.3 ) 
BACKPA8* F13,3 / 
BACKPAS* F13,3 ) 
PUZZLE * F13.3 / 
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• 6X * CP FOR PUZZLE • F13,3 ) 

600 FORMAT f 6X *NW8 * 110 f 6X *NWR * U0 ) 

700 FORMAT ( 6X *NBS * Xl0 / 6X *NUP * 110 ) 

cwn 

SUBROUTINE SETOATA 

common NOOF(2500)rIQ(20)fHE»PE(60)rST(60,60) 

COMMON / 0000 / MAXF,MAXB»MAXW,XPRNT 

COMMON / SUBT / N0EL*N0PT.XSAV,NSTR,MSU8,NUEL,NSU8,ISUB(63) 

COMMON / PRNT / XP0 # XP 1 , IP2, XP A , XPB, NUM I , NUM2 , IPS 

COMMON / MESH / IGO,NOFRE, NODES, NDXVY,NOXVX,NOXVZ 

COMMON / SCRl t IR,LP»Xl,X2,NR,MD,NW,S(204i) 

COMMON / TRAC / NHS,NWR,NBS,NUP 

LP s 0 

NWS S 0 

NWR B 0 

NBS 8 0 

NUP = 0 


DO 1 X 3 1,6 

1 CALL WIND ( X ) 

READ e00,MAXF,MAX8,MAXW 
PRINT R00,MAXF 
PRINT 901,MAXB 
PRINT 902,MAXW 

READ 100,IGO,XP0, XPl,XP2,IPA,IPB,NUMl,NUM2,IP8 
PRINT 100,XCO,XP0,IP1, XP2,XPA,IP8,NUMl,NUM2,IPS 
read 100,NOPT,NOEL»NDFRE,NODES,NDIVY,NOXVX,NOXVZ,IPRNT 

IP ( XG0,EQ,1 ) CALL SETUl 
IF ( 160, EQ. 2 ) CALL SETU2 
IF ( XG0,EQ,3 ) CALL SETU3 
IF ( XGO,EO,a ) CALL SETUfl 
IF ( IGO.EO.S ) CALL SETU5 

PRINT 200,NOPT ,NOEL »NDFRE 
PRINT 201, NODES, NDXVY,NDIVX 
PRINT 202,NDIVZ,IPRNT,LP 

CALL MIND ( 1 ) 

CALL WIND ( 2 ) 

100 format (20X4) 

200 FORMAT ( *1* 5X *NOPT * 15 / 6X *NOEL * 15 / 6X *NDFRE* 15 ) 

201 FORMAT ( 6X *N0DE3* 15 / 6X 8NDIVY* IS / 6X *NOIVX* 15 ) 

202 FORMAT ( 6X *NOIVZ* 15 / 6X *IPRNT* 15 / 6X ♦IPOS * 15 ) 

600 FORMAT ( 416 ) 

900 FORMAT ( *1* 5X *MAXF SET AT* 16 ) 

901 FORMAT ( 6X *MAXB SET AT* IS ) 

902 format ( 6X *HAXW SET AT* 16 ) 

RETURN 

END 

subroutine SETUl 

common NDOF(2500)f IQ(20)|HE,PE(60),ST(60,60) 

COMMON / 0000 / maxf,maxb,maxw,iprnt 

COMMON / SUBT / NOEL»NOPT, I8AV,N8TR,MSUB,NUEL»N8UB, ISUB(63) 
common / ME8H / IG0,NDFRE, NODES, NOIVY, N0IVX,NDIVZ 

READ 100,(NOOF(I),I81,NOPT) 
mo 


IF(IPRNT,GE,2) PRIMT IRQ# (WOOF(n,lBl,NOPT) 

DO 1 Ial,NOEL 

READ IM0.NODE8,ME, (10 (L) ,L* UNDOES) 
IFnPRNT,GE,2) PRINT 300 , NODES , HE .( I Q (L ), L = l» NODES ) 

CALL UNITl (l,t»NOOES,IQ) 

1 continue 

DO « L2 s UNOEL 

READ 100»ME 
IF(IPRNT,GE,2) print 100, me 

READ 500, (PE (M),HsI,mE) 
IF(IPRNT,GE,2) print 500, (PE (M),Ms1,mE) 

DO 2 X s t, ME 

READ 500,(ST(I,J),J=),ME) 

2 IF(IPRNT,GE,2) print 500, (ST( I , J) , J=1 ,ME) 

CALL UNIT2 ( ME, PE, ST ) 

U CONTINUE 


100 format (20141) 

300 FORMAT (6X,*JQ*,9I5) 

500 FORMAT ( 20F6.0 ) 

RETURN 

END 

SUBROUTINE 5ETU2 

COMMON NOOF(2500) ,10(20) ,ME,PE(60) ,ST(60,60) 

COMMON / SU8T / N(1EL,N0PT, ISAV,NSTR,MSUB,NUEL,NSUB,ISUB(63) 
COMMON / MESH / IG0,N0FRE, NODES, N0IVV,NDXVX,NDIVZ 

DO 1 I a l,NOP’f 

1 NOOF(I) a NOFRE 

DO 2 KalfNOEL 

read 130, nodes, me, (IQ(L) ,L= l ,NOOES) 

CALL READV 

CALL UNIT2 ( ME, PE, ST ) 

2 CONTINUE 

100 FORMAT (10I4J) 

RETURN 
END 

SUBROUTINE SETU3 
COMMON 

COMMON / 8UBT / 

COMMON / MESH / 

DIMENSION J0(4t) 


NDOF(2500),IQ(20),ME,PE(60),8T(60,60) 
NOEL»NOPT, I8AV,NSTR,M8UB,NUEL,N8UB,XSU8(63) 
IG0,NDFRE,N00E8,NDIVY,N0IVX,NDXVZ 


NOEL s NOIVY * NDIVX 
NOPT a (NOIVY+l)*(NOIVXfl) 

DO 1 X a 1,N0PT 
1 NOOF(I) a NOFRE 

ME a NDFRE * NODES 
NPl a NDIVY ♦ 1 
DO 3 IX a 1, NDIVX 

JOd) a (lX-1) * NPl ♦ I 
J0(41) s JQ(l) ♦ 1 
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PAGS 


IS 


Q^ALrry 


2 


JQ(2) » JQU) 4 NPl 

JGf3) s JQ(4) ♦ NPl 

00 3 lY s 1, NOIVY 

00 2 I a 1,4 

IQ(I) a J0(I) ♦ lY • 1 
CALL READY 

CALL UNIT2 ( ME»PE,ST ) 

3 CONTINUE 

RETURN 
ENO 

SUBROUTINE SETU« 

common NOOF (2500), 10(20), me, PE (60), ST (60, 60) 

COMMON / SUBT / NOEL,NOPT, ISAV,NSTR,MSUB,NUEL,N5UB, ISUB(63) 
common / MESH / XGO,NOFRE, NODES, NOIVY, ND1VX,N0IVZ 
OIMENSION JQ (6) 

NOEL « NOIVY * NOIVX 

NOPT a 3*NDIVY*N0IVX ♦ 2* (NDI VY+NDIVX) ♦ 1 

00 1 I a 1,N0PT 

1 NOOF(I) a NOFRE 

HE a NOFRE * NODES 
NPl B NDIVY ♦ 1 

INC ■ 3 * NOXVY ♦ 2 

JQ(1) B 1 

JQ(8) B 2 

JQCa) 8 3 

J0(5) 8 2 * NPl 

JQ(7) 8 2 * NPl t 1 

JQ(2) 8 3 * NPl 

JQ(6) s JQ(2) 4 1 

J0(3) 8 JQ(2) ♦ 2 

00 5 IX s 1, NOIVX 
LI a (9 
L2 so 

00 3 lY a 1, NDIVY 

00 2 I 8 1,4 

2 IQ(I) a JQ(I) ♦ L2 

10(5) a J0(5) ♦ LI 
10(6) a JQ(6) ♦ L2 
10(7) 8 JQ(7) ♦ LI 

IQ(S) 8 J0(6) 4 L2 

LI a LI 41 

12 s L2 4 2 

CALL ready 

CALL UNIT2 ( MEfREfST ) 

5 continue 

DO (I X 8 1, NODES 

4 JQ(X) 8 JO(I) 4 INC 

5 CONTINUE 

RETURN 
ENO 

SUBROUTINE SETU5 
COMMON 

common / subt / 

COMMON / MESH / 


NOOF (2500), IQ (20), ME, PE (60), ST (60, 60) 

NOEL, NOPT, ISA V,NSTR,MSUB, NUEL,NSUB, ISUB (63) 
IGO, NOFRE, NODES, NOI VY, NOIVX, NOXVZ 


dimension JQ(20) 

NOEL a MOIVX * NOIVV * NOTVZ 

Nl a (2*N0IVV4l)*(N0IVX*U ♦ (NOIVV^I) *ND1VX 

N2 3 ( NOIVYfl )*(NOlVX4l) 

N1N2 8 Nl 4 N2 
NnPT 3 Nl 4 NIN2 • NOIVZ 
DO I I 3 l,NOPT 

1 NOOF(I) 3 ndfre 

ME 3 NOFRE * NODES 
DO 2 I 3 1,3 

2 JQ( 1) a I 
1=24 34NOIVV 

DO 3 I 3 1,3 

3 JQCI43) 3 I 4 L 

DO <1 I 3 1,6 

4 JQ(l46) a JQ(I) 4 NIN2 

JO(l3) a 24NDIVY 4 2 
JO(l«) a JO(l3) 4 1 
JQ(15) a JQ(13) 4 N1N2 
JQ(16) 3 JQ(15) 4 1 
JQ(17) a Nl 4 I 

JQ(ia) a Nl 4 2 

JQ(19) a JQ(18) 4 NOIVY 
JQ(20) a J0(19) 4 1 

DO It IZ a I, NOIVZ 

U a 0 

12 a 0 

L3 a 0 

DO 9 IX a 1,NDIVX 

DO e lY a 1, NOIVY 

DO 5 I a 1,12 

5 lQ(t) a JQ(I) 4 Li 

DO 6 I a 13,16 

6 10(1) a JQ(I) 4 L2 

DO 7 I a 17,20 

7 IQ(I) 3 JQ(I) ♦ L3 

LI a tl 4 2 

L2 a L2 t t 

L5 a t3 ♦ I 

CALL READY 

CALL UNIT2 ( ME, PE, ST ) 

6 CONTINUE 

LI a ( 3*NDIVY42 ) * IX 
L2 a LI 
L5 a L3 ♦ I 
9 CONTINUE 

no 10 I a 1,20 

10 JOCn a JQ(I) 4 N1N2 

11 CONTINUE 
RETURN 
END 

SUBROUTINE READY 

common NOOF(2500),IQ(20),ME,PE(60),8T(60,60) 

COMMON / MESH / ICO, NDFRE, NODES, NOIVY, NDIVX, NOIVZ 
COMMON / 0000 t MAXF,MAXB,MAXW,IPRNT 

IF(IPRNT,GE,2) PRINT 100, NODES, HE 
IF(IPRNT,GE,2) PRINT 200, ( 10 (L) ,L«1 »NOOES) 
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CALL UNITl (!, 1, NODES, 10) 


00 I I 3 |,mE 

1 PE(I)s o.ia 

13 3 NOOF(l) 

L a e 

DO 2 I 3 I, ME, 13 
L a L ♦ 1 

K2 a I ♦ 13 - I 

00 2 K 3 i,Ka 

2 PE(K) 3 IQ(L) 


00 3 I a I, ME 
00 3 J a i,mE 
3 ST(I,J)3 1,0E*20 

00 tt I B i,ME 
« ST(I,X)a PE(I) 


100 FORMAT ( 6X 215 ) 

200 FORMAT ( 6X 1215 ) 

RETURN 

END 

SUBROUTINE UNIT2 ( ME, PE, ST ) 
dimension PE(60),STffe0,60) 



common 

/ 

SCRl / IR,LP,n,I2,NR, 


COMMON 

/ 

0000 / maxf,maxb,maxw. 


COMMON 

/ 

TRAC / NWS,NWR,NBS,NUP 


MO 

s 

ME 


NW 

s 

(M0*(MD^l))/2 ♦ MD ♦ 7 


IF 

( 

Nk»,GT,MAXB ) 


11 

s 

t 


12 

r 

MD 

1 

12 

s 

0 


IT 

a 

MO ♦ 1 

2 

11 

3 

12 ♦ 1 


12 

s 

n 


NW 

3 

7 


DO 3 I 

3 

n , MD 


NW 

IF 

3 

NW ♦ IT 
( nw.gt.maxr ) 


12 

3 

12 ♦ 1 

3 

IT 

a 

IT - 1 

a 

NW 

3 

NW • IT 

5 

12 

e 

12 - 1 

6 

NR 

3 

12 - 11 ♦ 1 


L 

8 

0 


DO 7 I 

3 

11 , 12 


L 

3 

L ♦ 1 

7 

S(L) 

3 

PE(I) 


00 8 I 

S 

11 r 12 


DO 8 J 

8 

I , MD 


L 

m 

L ♦ 1 

8 

S(L) 

3 

ST (I,j) 


IF( IS.NE.NU ) 


IR 

IS 

IR 


NW 

IR 

IR 


/ 6a 

* 6a 

♦ 1 

lOR 




GOTO I 
GOTO 6 

GOTO a 
GOTO 5 


NW 3 IR *64 • 7 

IP a IR ♦ LP 


CALL lOP (2HWB,P, IR, 7) 
CALL lOP (2HhB,2, S, KW) 
NWS 8 NWS ♦ NW ♦ 7 


IF ( IZ.LT.MD ) 


GOTO 2 


RETURN 

END 


subroutine FRONTIO 

COMMON 
COMMON 
COMMON 
COMMON 

common 
common , > 

data L1»L5 



LQK950) 

lfNUEL>NSUB»I8UB(63) 


READ S,MSUB 
PRINT I0,M8U0 
IM2 3 MSUB ♦ 1 


CALL SECOND ( T1 ) 


ISAV 8 NOEL 
NSTR a 0 
DO 1 IM 3 1,IM2 
NT 3 0 

1 call SUBIOS ( NT,IM ) 
DO 2 IM 8 1,6 

2 CALL HIND ( IM ) 


5 FORMAT ( la ) 

IP PORMAT ( *1* 5y 

return 

END 

SUBROUTINE 
common 

COMMON 

common 

common / subt 
common / 0000 


♦NUMBER OF SUBSTRUCTURES IS* 15 ) 


data Ll.LS /I, 5/ 


SUBIOS ( NT,XM ) 

NOOF (2500), IB (2500), IC (2500), LO( 150), LOC( 150) 
NODEl, 101(150), Ml, Kl,Nl,NOl,MOCl (150) 

NOOE2r IQ2( 150) ,M2#K2»N2f N02#MOC2( 150) 

NOEL# NOP T # ISAV,NSTR,MSUB,NUEL#NSUB, ISUB(63) 
maxf,maxb,maxw,iprnt 


/ 

/ 


IF (IM,LE,MSUB) 

CALL LOCATE ( NOEL # IM, NT,L1 ,L 1 ) 


GOTO 1 
GOTO 5 


t READ 

IF(IPRNT,CC,2) PRINT 

CALL NRIT3 (1) 


10i NSUB#(I8UB(I),l3l,NSUB) 
20 #XM,NSUB« (ISUB(I),Is1,nSUB) 


DO a 


NUEL a 0 
XT 8 l.NSUR 

IF(IPRNT*GE,2) 


READ 10#NODE2, (I02(L) #L8l ,N qdE 2) 
PRINT 30,NODE2, (I02(L) >Lal ,N00E2) 


00 2 I 8 1,N00E2 
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, L = IQ2(I) 

2 NOOF(L) s»nOOF(L) 


ISUB a TSUB flT) 

CALL WIND ( 6 ) 

call locate t LSU0,IM,NT,Lt»6 ) 
NUEL B NUEL ♦ LSUR 


00 3 la 1 ,nopT 

NOOFU) 3 ZABSCNOOFCU) 


CALL UNIT! ( L5#1,N2»M0C2 ) 

CALL MERGEIO ( L t »L5, N0DE2, IQ2) 


10 FORMAT ( 2Pia ) 

20 format (*0* *IM 
30 FORMAT ( 6X*M0 


*22Ifl ) lU / 6X AISUB* 2014 ) 


i RETURN 
END 


COMMON / SU8T / NOEL,NOPT,I8»V,NSTR,MSUB,NUEL.NSUB,ISUB(6J) 


NOEL = NOEL - NUEL 
IF t NOEL.EQ.B ) 


00 1 X a 1 ,noEL 


CALL UNITl ( LI »2, NODES, lO ) 
call UNITl ( L5, I, NODES, 10 ) 

1 CONTINUE 

2 NOEL a NOEL ♦ NSUB 

NSTR a NSTR ♦ NSUfl 

call wind ( LI ) 

CALL WIND ( L5 ) 

LT a LI 
11 « L5 

L5 a lT 


return 

end 


subroutine 
common 
common 
common 
common 


/ 8U0T / 


locate ( L8UB,IM,NT,LI,LT ) 

NOOF (2500) , ZB (2580) , IC (2500 ) ,L0( 150) «LOC( 1 RBI 
NOOE1,IO|(,5«),M1,K|,n,.noWMOC1MSB) ® 
node 2# 102(150) ,M2,K2,N2,NQ2,H0C2( 150) 
NOEL,NOPT,I8AV,N8TR,M8U0,NUEL,NSUB,I8UB(63) 


NQL a 0 


00 1 Z ■ l,N0PT 
1 1B(1) a 0 

00 R LZ ■ 1,LSU8 


original PAoe IS 

OF ntmitx*. 


CALL UNIT! ( Ll,2,NO0E2,lO2 ) 
XF (IM,CT,MSUB) 
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GOTO 2 


C4LL UNITl ( LT»l,NO0a2,I02 ) 

2 00 3 1 a 1,N00E2 

K » 102(1) 

3 IB(K) a I8(K) ♦ 1 

« CONTINUE 


00 5 X a |,N0PT 
5 XC(I) a XB(X) 


CALL WIND ( LT ) 

LB a LSUe • 1 

CALL front ( NT,LT,0, 0 

IF ( LB,EO.0 ) 

00 6 LZ a 1,LB 

KZ 

CALL UPDATE 

CALL front ( NT,LT,0,KZ 

6 CALL EXPAND ( LZ 

7 call update 

CALL FRONT ( NT,LT,1, 0 
CALL EXPAND ( LSUa 


) 

GOTO 7 

a LB « LZ 

) 

) 

) 

) 


RETURN 

END 


subroutine front 

common 

common 

common 

common / SUBT / 


( NT,LT,NG,KZ ) 

NOOF (2500), IB (2500), XC (2500), L0( 150), LOC 050) 
NC0E1,IQ1 O50),M1,KJ,NI,NQ1,MOCI (150) 

N0DE2, 102050) , M2, K2,N2»N02, M0C2( 150) 
NOEL,NOPT , ISAV,NSTR,MSUB,NL'EL,NSU5, I8UB(63) 


IF ( NG.GT.0 ) 

CALL UNITl ( LT,2,N00E2, IQ2 ) 


GOTO 2 


DO 1 I a 1,N00E2 

K a IQ2(I) 

IF ( NDOF(K),GT,0 ) 
IF ( NDOF(K),LT,0 ) 

I continue 


IC(K) a IC(K) - 1 
IC(K) a IC(K) ♦ 1 


2 n2* 0 

00 3 I a l,NOPT 

IF ( 10(1), EO, ICO) ) 

IF ( ICO), EG, 0 ) 1H(I)80 

W2 a n2 ♦ I 
L0C(N2) a I 

3 continue 


GOTO 3 


IF ( N2,EQ,0 ) 

H2 a 0 

00 4 I a i,n 2 
MOC2(I)a0 

J ■ LOCO) 

IF ( IC(J),NE.0 ) 
LOC(I)80 
M2 a M2 ♦ 1 
M0C2(M2) a j 
4 CONTINUE 


GOTO 6 


GOTO 4 
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K2 s H2 ♦ 1 

CAUL gaps ( NT,KZ ) 


6 RETURN 
END 

subroutine gaps ( NT»KZ ) 

CO^'BON NOOF (2500), I0(25B0),IC(25CP),LO(l5j!)),LOC (150) 

COMMON MOOE1#XOI(150) ,Ml,Kl,Nl,NOl,MOCl (150) 

COMMON NODE2,IO2(150),M2,K2,N2,NQ2,MOC2nS0) 

COMMON / 0000 / MAXF,MAXB,HAXW,XPRNT 


IF ( KZ ,EQ,0 ) GOTO 6 

IF ( NQl.GT.MAXF ) GOTO 0 

IF ( NQ2,GT,MAXF ) GOTO 0 

1 It s HAX0 ( K1,K2 ) 

12 B HIN0 ( N1,N2 ) 

IF (II, CT, 12) GOTO 10 


00 U I3Z1,I2 

KbMOCKX) 


00 2 J31,N2 

IF ( LOC(J),EQ,K ) goto 3 

2 CONTINUE 


3 

4 


M0C2(I)8K 

LOC(J)B0 

continue 


uOTO 4 


DO 7 Ial,N2 

IF ( LOC(I),EO,0 ) GOTO 7 

CO 5 JBK2,N2 

IF ( MOC2(J),EQ,0 ) GOTO 6 

5 continue 

6 M0C2(J) 8 LOC(I) 

UOC (1)30 

7 CONTINUE 

GOTO 10 

0 J8K2 

00 9 X«l,N2 

IF ( LOC(I),EQ,0 ) GOTO 9 

M0c2(J)3L0C(I) 

LOC(I)30 

J » J ♦ I 
9 CONTINUE 


10 NQ2 I NUM ( N2,M0C2,N00F ) 

NT B HAX0 ( NT,NQ2 ) 


RETURN 

END 

subroutine UPDATE 

COMMON NOOF (2500) , IB (2500) , XC (2500) ,LQ( 150) ,LOC( 150) 

common NOOElflQl (150) ,M1 ,k1,N1,NQ1,H0C1 (150) 

COMMON N0DE2, 102(150) , M2, K2,N2,NQ2,M0C2( 150) 

00 1 I 8 1,N0DE2 

1 IQl (I) 8 102 (I) 
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00 3 Is 1,N2 

3 MOCl(t) a M0C2(I) 


NOOEl 

8 

N0DE2 

Ml 

a 

M2 

K1 

8 

K2 

Nt 

8 

N2 

NQl 

8 

NQ2 


RETURN 

END 

subroutine expand ( LZ ) 

common MOOF(25eO),TB(2500),IC(250D),LO(l50)»LOC(l50) 

common NOOEWIQI (150),MI,KI,N1,NQ1,MOC1 (150) 

COMMON MOOE2»ZO2(150) »M2#K2,N2,NO2#MOC2(150) 

COMMON / FRNT / LO,MO,KO,NQL,NQ,NOR,MOI»LOT(950) 

DO 1 I 8 IfNOOEl 
00 I J a l,Nl 

1 IF ( XQKn.EO.MOCKJ) ) LQ(I) 8 J 
LO 8 0 

00 2 I 8 l,NOOEl 
tP » lOl ( I ) 

J2 8 NUM ( LO(I)>HOCtf NDOF ) 

J1 8 J2 - IABSCNOOFCUP) ) t 1 
DO 2 J 8 Jl,J2 
CD 8 LO ♦ 1 

2 LOT<L0)aJ 

MOl 8 NUM ( N0DE1|I01,ND0F ) 

MO e NUM ( Ml,MOCl,NOOF ) 

KQ 8 MQ41 
NO 3 NQt 

IF ( Kl.GT.Nl ) goto B 

00 5 I 8 Kt»Nl 
00 3 J 8 W N2 

IF ( M0C1(X),EQ.M0C2(J) ) goto Q 

3 continue 

LOC ( I ) 8 J 

5 CONTINUE 

IF ( MQ,Eo,0 ) GOTO 9 

DO 6 I 8 i,MO 

LD 8 to ♦ I 

6 LOT(tO) 8 0 

9 00 7 I*Kt#Nl 

LP 8 MOCl ( I ) 

J2 8 nUm ( L0C(I)«M0C2rND0F ) 

■ J2 • IABS(N00F(LP)) ♦ 1 
DO 7 J8J1,J2 
L0«L04t 

7 L0T(LD)8j 


S 


NOR 8 N02 


ni 


CALL WRIT3 (2) 
NOL 5 NQ 


CALL PRNTiJ ( LZ ) 

CALL PRNTt 
call PVNT2 ( LZ ) 

RETURN 

END 

FUNCTION NUH ( lPOS»HOC»NDOF ) 

dimension Mocn)*NDOFn) 


NUM s 0 

IF ( IPOS.EO.e ) 

DO 1 I a J,IPOS 

K 8 HOC ( I ) 

1 NUM a NUM ♦ lAflS ( NOCF ( K 

2 RETURN 
END 

subroutine FRONTS! 

COMMON 

COMMON / FRNT t 
common / SUBT / 
common / 0000 / 

COMMON / WRIT / 


GOTO 2 


) ) 


8(700), Adseee) 

LO,MO,Kn,NQL ,NO,NQR,MOI ,LOT (950) 
NOEI.,NOPT,ISAV,NSTR,MSUB,NUEL,NSU£ 
MAXF,MAXB,MAXW,IPHNT 
LPU9,NBAC,Nwi,nw2,NC69,N0L0 



LPU« 

a 0 


NOLO 

a 0 


1M2 

a msuB ♦ 


NOEL 

a ISAV 

DO 3 

IM a 

1,IM2 


CALL 

SUBSTR ( 

DO « 

IM s 



IF ( 

XH.NE.O ) 


( NT,IM ) 


RETURN 
END 

subroutine SUP.STR ( NT,IM ) 
COMMON 

COMMON f fRNT / 

COMMON / SUBT / 


B(700),A(U000) 

LO,M(3,ko,NOL,NQ,nOR,MDI ,LOT (950) 
^0EL,N0PT, lSAV,NSTR,MSUB,NUEL#NSUe 
COMMON / 0000 / MAXF,MAXB,HAXW, IPRNT 
data 


IPOS a 0 

DO 1 I a 1, 700 

1 B(I) a 0 

DO 2 I a 1,13000 

2 A(I) a 0 

IF ( IM.LE.MSUB ) 

CALL 8CANI0 ( NOEL, LUL2, IPOS ) 

3 CALL WIND ( LI ) 


GOTO 3 
GOTO 5 


,ISUB(63) 


,ISUB(63) 
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call REA03 (1) 


( 


NUEL a ft 

DO a IT = i,NSU8 

LSUfl a I8UB ( IT ) 

f;*LL SCANIQ ( LSUB,L1,L2,IP0S ) 

NUEL 3 nuEL ♦ LSU8 

CONTINUE^^'”^^^ ^ ^ 


U*1P0S ) 


CALL WRIT« M,NE) 


CALL MERGEST ( U,L2#IP0S ) 
IF(IPRNT,GE,2) print 20,IPO3 


20 


format ( 6X *IPOS IS * 14 ) 
return 

END 


SUBROUTINE MERGEST ( LJ,L2,IP08 ) 

io.z ; s*c“p^i ; 


NOEL s NOEL •• NUEL 
IF ( NOEL,EO.0) 

00 2 I s I, NOEL 


GOTO 3 


» CALL SEMBl ( L2»2,IP0S ) 

CALL SEM01 ( Ll»l»IPOS ) 

IF ( 12,LT,H0 ) 

2 CONTINUE 

3 NOEL s NOEL ♦ NSUB 

call NINO ( LI ) 

CALL WIND ( L2 ) 

LT a Lt 
LI s L2 
L2 s lT 


GOTO I 


return 

END 

‘ t-Sueai.LZ.IPOS ) 
common B(700),a(130B0) 

COMMON / FRNT / LO.HQ . KQ . Nm . kin _ kino I 


w V r • ^ 

common 

COMMON / 

common / 

COMMON / 
INTEGER 
common / 
LOGICAL 
DATA L5 


• — ^ ^ r w * f I. f A r u w / 

B(700),a(13000) 

LO,MQ,KQ,nQL»NQ,NQR,M01 ,LQTC950) 
R1»R2,RM,WM,H1,H2,NE#NC 
MAXF,MAX8,MAXW,IPRNT 
R1»R2,RM,WM,H|,H2 
!22 / LOO,LOLO(700) 

Lo / 5 » 6 / 


FRNT 

1111 

0000 

2222 
TE 
# L6 


00 2 LZ • 1,L8UB 


UDU 


= [a 


CALL REA03 (2) 

CALL FORMNO ( NO ) 


TE8TLF a (NQL.GT.MAXF ^,0R, NO ,GT,MAXF) .AND, LZ.CT.l 


2 ^ 


■ I I 


TESTRF 

s 

(NQ .GT.MAXF 

.OR, 

NUR, 

GT, 

HAXF) 

DI8KLT 

3 

TESTRF 

.AND, 

LZ 

.LT 

, LSUB 

DISKEQ 

3 

TESTRF 

.AND, 

LZ 

• EO 

, LSUB 

SETCOR 

S 

.NOT, TESTRF 

.and. 

LZ 

.LT 

. LSUB 



JPOS 

a 0 






NB 

a 0 






NR 

a NB 

♦ 1 





CALL 

ROMS 

C 

NB 

) 


CALL PPNTSn,LZ>NB, NOTH, NOTH, LP) 


CALL SEMBLEI ( NB,L2 ) 


IF 

( 

TESTLF ) 

CALL 

SEHBLE2 

( NB,L5 ) 




CALL 

PRNTA 

(A,6,R2,NQ) 

IF 

( 

NB,EQ,1 

) 

Rl a 2 


IF 

( 

HQ,GT,0) 

CALL REDUCE ( 

nb,mq,nq ) 

IF 

( 

NB,EQ,1 

) 

Rl a KQ 


IF 

( 

OISKLT ) CAUL 

SETDISK 

( L6,JP0S ) 

IF 

( 

OISKEQ ) CALL 

SETOISK 

( Ll.IPOS ) 

IF 

( 

R2,LT,NQ ) 


GOTO 1 




CALL 

WRITfl 

(2.NE) 




CALL 

ZEROA 

(NE,L5,L6, TESTLF, TESTRF ) 

IF 

( 

SETCOR ] 

I CALL 

SETCORE 

( NQ,NQR ) 

IF 

( 

TESTRF ] 

1 CALL 

PRNTS 

( 2»lz»noth,jpos»noth,lp ) 

CONTINUE 





RETURN 

ENO 

subroutine ZEROA (NE,L5,L6»TESTLF»TESTRF) 

COMMON B(7fl0),An300U) 

COMMON / FRNT / LD,MQ,KQ,noL,NQ,N(3R,MOI,LQT(950) 
COMMON / 2222 / LOO,LOLD (700) 
logical TESTLF,TE8TRF 

IF ( MQ,EO,0 ) goto 3 

DO 1 I 8 1,MQ 

1 BCD s 0 

DO 2 I s l,NE 

2 A(I) 8 0 

3 LOO » 0 

IF (.NOT.TESTRF) GOTO 5 

L a 0 

DO a I s KO,NQ 
L a L ♦ I 

tt LOLD (L) a LOT (I+MDl) 

LOO a L 


LT a L5 


n ^ 


IS 


L5 

L 6 


L 6 

LT 


5 IF ( TESTUF.OR.TtSTRF ) CALL WIND ( L5 ) 

IF ( TESTLF.OR.TESTRF ) CALL WIND ( L6 ) 


RETURN 

END 

SUBROUTINE ROWS 
COMMON / NOOO / 
common / FRNT / 
common / 00B0 / 

COMMON / nil / 
INTEGER 

COMMON / 2222 / 


( NB ) 

NO(700) 

LD,MQ,KQ,NOL,NO,NO«,HOULOT(950) 
MAX^,MAX0,MA)(W,IPRNT 
R1 f R2 ,Rm,wm,H1 ,H2#NE,NC 
R1 #R2»RM, WM,H1 ,H2 
LOO,LOLO(700) 


J 

2 


DO 


3 

5 


3 


IF ( NB.GT.l ) goto 1 

Rl = I 

R2 = MO ♦ 1 

RM a 0 

WM s 0 

NE 8 0 

IF ( MQ,GT,0 ) NE s NO(MQ) ♦ NO - MO 

IF ( NO(NO),GT.HAXW ) GOTO 2 

R2 s NO 

NC = NO(NO) - NE 


Rl 8 R2 ♦ I 
R2 8 Rl 
NC 8 NE 

IT 8 NQ - MQ - RM 
II = R2 
I a II, NO 
NC 8 NC ♦ IT 
IF ( NC.GT.MAXW ) 
R2 8 R2 t I 
IT = IT ip I 

NC a NC - IT 
R2 8 R2 • I 


GOTO 6 


GOTO 4 
GOTO 5 


NC 8 NC «* NE 

6 HI 8 0 
H2 a 0 

DO 8 I 8 R1,R2 
DO 7 J 8 I, MO I 

IF ( LOT (J),NE,I 
HI 8 HI ♦ 1 

7 CONTINUE 

6 continue 


GOTO 7 
GOTO 8 


DO 10 
DO 9 


IF ( LDO,EO,0 ) 

I 8 R1,R2 
J a I, LOO 
IF ( LOLD(J),NE,I ) 
H2 8 H2 ♦ I 
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GOTO 11 


GOTO 9 


CONTINUE 

CONTINUE 


GOTO 10 
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U RETURN 
END 

SUBROUTINE SEMBLEl ( NB,t2 ) 

COMMON B(7BU) f A(130BR) 

COMMON / NOOO / NO(70O) 

COMMON / FRNT / LD,MQ,KO»mqu,NQ,NQR,MDI,LOT(950) 

COMMON / SCRl / IR,LP#Il»I2.NR,MO#NW,S(20«l) 

COMMON / tm / Rl,R2,RM,WM,Hl,H2,NE»NC 

INTEGER Rl,R2*RM,wM,Ht»H2»R ,BACK,FINI 

IF ( HI ,EO.0 ) GOTO 5 


IF ( NB ,EO,l ) LOGO s 0 

LOGR = 0 
BACK = 0 
FINI a 0 
I dec' a 0 

1 IF (LOGO.NE.n CALL SEMBl (L 2 , 2 ,N 0 TH) 
LOGO a 0 

CALL PRNTS(3»N0TH,N0TH,N0TH,I0EC»LP) 

R a 0 

t a NR 

DO 3 M a 11,12 
R a R ♦ 1 
I a LOT(M) 

K a LQT(M) • RM 
IF ( I,GT,R2 ) 

IF ( l.LT.Rl ) 

IF ( 1,6T,R2 ) 


2 DO 3 N a M , MO 

I a MIN0(LQT(M),LOT(N)) 

J a MAX0(LQT(M) ,LOT(Nn 
K a NO(I) • I ♦ J - HM 
L a L ♦ I 

IF ( I.LT.Rl ) GOTO 3 

IF ( I.GT.R2 ) GOTO 3 

A(K) a A(K) f S(L) 

3 Continue 

IF ( BACK<fFINI,EO.0 ) GOTO 4 

LOGR a LOGR 4 1 
RACK a back * IR 

a IF ( HI ,GT,0 ) GOTO I 

CALL PRNTS (<»» noth# NOTH, NOTH, NOTH, LP) 

LOGO a LOGR 

IF (LOGO.GT.l) CALL SEMBl (L2,3,BACK) 

5 RETURN 
END 

subroutine SEMBLE2 ( NB,LS ) 

COMMON B(700)#A(13000) 

_ 1 1 ^ 


FINI a 1 

GOTO 2 
GOTO 2 

Ml a HI - 1 
B(K) a b(K) t SCR) 


COMMON / NOOO / 
COMMON / 2222 / 
COMMON / SCR2 / 
common / nil / 
INTEGER 


NO(70t!) 

LOO»LOUOC700) 

IR»LP» 1 1 » I2,NR,MD#NW,S(20ai ) 
R1,R2 ,Rm,mm,hi,H2,NE,NC 
R1,R2,RM,wm,HI,H2,R »BACK,FINI 


IF ( H2 *£0,0 ) GOTO 5 


IF ( NB (EOtl ) LOGO B 0 

LOGR e 0 
BACK a 0 
FINI = 0 
XOEC a 0 


I IF (LOGO.NE.l) CALL 5EMB2 (L5,2,N0TH) 
LOGO a 0 


CALL PRMTS (3, NOTH, NOTH, NOTH, IOEC,LP) 


R a 0 
L a NR 

DO 3 M a 11,12 
R a R ♦ 1 
I a LOLDCM) 

K a LOLn(M) * RM 
IF ( 1,GT,R2 ) 

IF ( I.LT.RI ) 

IF ( I,GT,R2 ) 


FINI a I 

GOTO 2 
GOTO 2 

H2 a H2 « 1 
B(K) a B(K) ♦ S(R) 


2 


3 


DC 3 N a M , MD 


I a min 0(LOLO(M),LOLOCN)) 
J a HAX0fLOLO(M),LOLO(N)) 
K a NOm • I ♦ J - MM 
L a L ♦ 1 
IF ( I,LT,R1 ) 

IF C I,GT,R2 ) 


CONTINUE 


A(K) a A(K) ♦ 8(L) 


GOTO 3 
GOTO 3 


U 


IF ( 8ACK+FINI,EO,0 
IF ( H2 ,GT,0 ) 


GOTO U 

LOGR a lOGR ♦ 1 
BACK a back f IR 

GOTO 1 


CALL PRNTSCa, NOTH, NOTH, NOTH, NOTH, LP) 


. * logo a LOGR 

IF (LOGO.GT.l) CALL 8EM02 fL5,3,0ACK) 

5 RETURN 
END 

subroutine reduce ( NB,MO,NQ ) 

common B(700),A(13000) 

common / NOOO / NOC700) 

common / nil / RI,R2,RM,WM,HI,H2,NE,NC 

INTEGER RI,R2,RM,WM,HI ,H2 


DO 1 I B R1,R2 


Lb I - RM 

XO « NQ(I) « I - MM 
J2 a I • I 
IP ( I.GT.Mq ) J2 3 mo 

DO 1 J a i,J2 

JO a NOrj) » J 
C B«A(J04I)/A(.I04J) 

B(L) B B(L) ♦ B(J) * C 

DO 1 K a X,NQ 

I A(|0>K) a A(XOfK) ♦ A(JO^K)*C 

RETURN 

END 

SUBROUTINE SETDISK ( L6#LP08 ) 

COMMON B(700)f AC13B00) 

COMMON / NOOO / NO(700) 

common / FRNT / LOrMQ,KO»NOl.,NO,NQR,MOULQT(950) 
common / STOR / IR,LP#JI#J2#NR,M0,NM 
COMMON / IIU / Rl,R2,RM,WM,Hl,H2,NE,NC 
common / 0000 / MAXF»MAXBrMAXW,IPRNT 
INTEGER RlfR2,RM,WM,Hl,H2 

MD a (go 1 MQ 


12 a Rl • 1 

IT a NO «• MQ » RM ♦ 1 

1 II a 12 ♦ 1 
12 a H 

NW a 7 

DO 2 X a It f R2 
NM a NW 4 IT 
IF (NM.GT.MAXB) 

12 a 12 4 1 

2 IT a IT . 1 

3 NM a NW • IT 

ft 12 a 12 • i 

NR a 12 • It ♦ I 


J1 a n • MQ 
J2 a 12 . MQ 
Kt a It i. RM 
K2 a 12 ^ PH 
Lt a NO(It) • WM 


L2 a NOCI2) ♦ NO • 12 • WM 
call store ( L6,LPOS#Kt,Lt,0»A ) 

DO 5 I a K1,K2 

5 8(1)30 

DO 6 I a LlfL2 

6 A(I)a0 

CALL PRNTS( 5, NOTH, NOTH, NOTH, NOTH, LP) 


GOTO 3 
GOTO 4 


IIB 


IF (I2,LT,H2) 


GOTO 1 


RM s RM ♦ R2 - Rl ♦ I 
MM S MM ♦ NC 

RETURN 
END 

subroutine SETCORE ( NO, NR ) 

COMMON e(7?D) r A(13B00) 

COMMON / NOOO / NO(7a0) 

IF ( NR.EQ.B ) 

IF ( NR • NO ) 

1 CALL SCRMBLF ( NO ) 

L s NR 

00 2 I s 2#NR 

N s NO(I) - I 

00 2 J s I, NR 

L 3 L ♦ I 
A(L) = A(N+J1 

2 A(NfJ) = 0 

CALL FORMNO ( NR, NO ) 

GOTO 7 

3 IN s NR - NQ 

L2 3 (NR*(NR + n )/2 - (IN*(lN + l))/2 ♦ I 

no 5 II 3 l,NQ 

I s NO - II ♦ 1 
IN a IN ♦ I 
LI 3 NOCn ♦ II 
L2 3 L2 ^ In ♦ II 
DO a L 3 1,11 

TEMP a A (Ll^L) 

A (LI-LT 3 0 
U A (L2-L) a TEMP 

5 L2 3 L2 • II 

CALL FORMNO ( NR, NO ) 

6 CALL SCRMBLE ( NR ) 

7 RETURN 
END 

subroutine SCRMBLE ( NT ) 

COMMON / FRNT / LO , MQ , KQ , NQL » NO, NOR , MO I ,LOT (950 ) 
KKOUNT 3 0 

IF (NQ ,GE, NT) GOTO 2 

J a NO 4> 1 

DO I I 3 J,NT 

1 lot (NMDl) a 0 

2 DO <1 M 3 KQ^NQ 

3 N 3 LQT(MfMOl) 

IF ( N ,EQ, M ,0R, N ,EQ, 0) 

. . nn 


GOTO 7 
I » ^ » 3 


GOTO 4 


LQT(M^mdI) 8 UOT(N+mDI) 
LQT(NiMOl) a K 


I s M1N0(M,N) 

J 8 HAXO(H,N) 

CALL SWITCH tI,J,NT) 


KKOUNT 8 KKOUNT f 1 
IF (KKOUNT ,GT, NT) STOP 

ti CONTINUE 


GOTO 3 


RETURN 

ENO 

subroutine switch ( M,N,NT ) 
COWHON B(7R0)#A(130t}0) 

COMMON / NOOO / NO(7VJ0) 


C =B(M) 

B(M)s8(N) 

B(N)sC 


HSsNO(H)>«H 

NS8NO(N)>N 


C 8A(MS^H) 
A(HSfH)3A(NSfN) 
A(N$4 N)sC 


IF ( M.EO.l ) GOTO 2 

I23M-1 

00 I 1= 1,12 

ISbnoU)-! 

C BAdS^M) 

A(IStM)sA(IS4N) 

1 A(IStN)8i: 


2 JlsM+i 

IF ( Jl.EQ.N ) GOTO 4 

J2sN-l 

DO 3 JBJ1,J2 

J8bN 0(J)«J 
C 8A(MSfJ) 

A(mS«J)8A(JS^N) 

3 A(JS^N)bC 


4 IF ( N.EQ.NT 
JtsN+l 

00 5 J8J1,NT 

C BA(MSfJ) 

A(MS^J)aA(NS4j) 

5 A(NS«J)8C 


GOTO 6 


OWGlNAl 


PAGE tS 


b RETURN 





END 

subroutine backsub 

COMMON H(7(^0) ,A(13B4»0) 

COMMON / SUHT / NOEL»NOPT,ISAV,N8TR,NSU0,NUEI»NSU8,ISUB(63) 
data L5 » L6 / 5 » 6 / 


CALL REAOO (3) 

IF ( HSUfl.GT.O ) 

call solve ( U l»NOEL»Lb ) 

1 left n NOEL - 1 

CALL SOLVE ( 1, i,LEFT,L6 ) 

CALL SOLVE ( 1* 2* 1»L6 ) 


GOTO 1 
GOTO U 


DO 3 IM a I^MSUB 


CALL WIND f L5 ) 
CALL WIND ( Lb ) 
LT a L5 
L5 a Lb 
Lb a LT 


CALL REA04 (1) 


DO 2 IT a IfNSUB 

LT a nSUB - it ♦ 1 
LSUB s ISU0 ( LT ) 

CALL UNITED (L5,2,MD,0) 

CALL PRNTB ( B,MO ) 

2 CALL SOLVE ( IM,MSUB»LSUB,Lb ) 

3 CONTINUE 


a RETURN 
END 

subroutine solve ( IM, JM,LSU0,Lb ) 
common B(7P0) ,An3000) 

COMMON / FRNT / LD,Mq,KO,N(3L,NQ,NQR,MOi ,LOT(950) 

COMMON / SCRl / IR,LP,Il,I2,NR,MD,NW,S(20tiU 

COMMON / SUbT / NQELf N0PT,ISAV,NSTR,MSU8,NUEL>NSUB,ISUe(b3) 

IF ( LSU0.EQ.0 ) GOTO 0 

DO 7 LZ a i,L8UB 

CALL READO (2) 
call PRNT2 ( LZ ) 

call FORMNO ( NQ ) 

IF ( KQ.GT.NQ ) GOTO 2 

DO 1 I « KQ,NO 

L a lot ( I ♦ MDl ) 

1 S ( I ) a B ( L ) 

2 DO 3 1 a 1,NQ 

3 B ( I ) a SCI) 

CAUL PRNTB ( B,NQ ) 
call PRNTA C A,B,MQ,N0 ) 

IF C MQ.GT. 0 ) call brass ( HQ,KQ,NQ ) 
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CALL PRNTB ( B,NQ ) 


« 00 5 I 8 1,M01 

^ L » LOT ( 1 ) 

5 8 ( I ) n B ( L ) 

IF ( ) 

CALL PRNTEO ( 2»ISAV,M0l,S ) 

CALL UNITED ( \, l^HOl^S ) 

6 CALL PRNTEO ( I#NSTR,MD1,S ) 

CALL UNITED (L6, t,H01,S ) 

7 CONTINUE 

8 RETURN 
END 

SUBROUTINE BRASS ( MQ,kO,NQ ) 

B(700),A(13BB0) 
COMMON / nOOL / NOf700) 

IF ( MQ ,E0, NO ) 

00 1 I s i,MQ 

K a NO(I)-I 
00 1 J a KO^NQ 

1 B(I) 8 B(I)«Af J4K)*B(J) 

?. Ha MQ4l 

I « NO(MQ) 

B(mq) 8 B(M0)/A(1) 

IF ( MO.EO.l ) 

00 « L s 2,MQ 

I a M-L 
JJ 8 !♦! 

K « NO(I)*I 

DO 3 J 8 J1,MQ 

3 B(I) a B(I)-A(JtK)*p(J) 

K 8 NO(I) 

« B(I) 1 B(I)/A(K) 

5 RETURN 
END 

subroutine FORMNO ( NT ) 

COMMON / NOOO / NOf700) 

Nocn 8 1 

L ■ NT+2 
00 I I 8 2, NT 

NO(I) 8 NO(I-l)fL - I 

return 

END 

subroutine UNITI ( NTAPE, ICO, NODES, 10 ) 
dimension 10(1) J 

17 


GOTO 6 
GOTO 7 

GOTO 2 
GOTO 5 


COMMOW f TRAC / NWS,NwR,NRS,NUP 

GOTO (1,2) IGO 

1 CAU lOP (2HWB,NTAPE, NODES. 1) 

CALL IDP (2HWR,NTAPE, IQ .NODES) 

Nwa s Nwy ♦ nodes ♦ I 

GOTO 3 

2 CALL lOP (2HRB,NTAPE, NODES, 1) 
call IOP (2HRB,NTAPE, IQ , NODES) 

NWR s NWR ♦ NODES t 1 

3 RETURN 
END 

subroutine WRIT3 ( IGO ) , • 

COMMON / TRAC / NWS,NmR,N3S,NUP 

COMMON / SUBT / N0EL.N0Pr,ISAV,NSTH,MSUB,NUEL»NSUB,I8UB(63) 
common / FRNT / LD,mq,kq,NQL,NQ,NQR,MOI ,LQT(R50) 

DIMENSION I0UM(956) 

equivalence ( MQ,I0UM ) 

GOTO ( 1,2 ) IGO 

1 CALL IOP (2HW0,3,NSU0,6a) 

NWS s NWS ♦ 

GOTO 3 

2 LO s LD ♦ 6 
CALL IOP (2Hwe,3,LO , r 
CALL IOP (2HWB,3,IDUM,LD) 

NWS B NWS 4 LD ♦ 1 

NUP 3 NUP 4 (HQ*(M0*MQ434NQ*NQ-3*MQ*NQ»l))/b 


3 RETURN 
END 

subroutine REA03 ( IGO ) 
common / TRAC / NWS, NWR, NBS, NUP 

COMMON / SUBT / N0EL»N0PT,ISAV,N8TR,MSUH,NUEL,NSUB,ISUB(63) 
COMMON / FRNT / LO,MQ,KO,NQL,NQ,NQR,MD1,LQT(950) 

DIMENSION IDUMC956) 

equivalence ( MQ,IOUM ) 

GOTO ( 1,2 ) IGO 

1 CALL IOP (2HRB,3,NSUB,69) 

NWR 3 NWR 4 64 

GOTO 3 

2 CALL IOP (2HRB,3,L0 , 1) 

CALL IOP (2HRp,3,10UM,L0) 

NWR s NWR 4 LO 4 1 


3 


I 


RETURN 

END 

SUBROUTINE WRIT4 ( IGO,NE ) 

COMMON B(7HM) » A(130P'J) 

COMMON / TRAC / NWS, NWR, NBS, NUP 
COMMON / FRNT / Ln,MQ,KO,NOL,NQ,NQR,MDl,LQT(950) 
common / SUBT / N0EL,N0PT,1SAV,NSTR,MSU6,NUEL,N8UB,ISUB(63) 
COMMON f WRIT / LPU4,NBAC,NW1 ,NW2,NC64,N0L0 
DIMENSION I0UM(956) 

equivalence ( M0,10UM ) 


GOTO ( 1.2 ) IGO 


IR =1 

LPUttsIR4LPU4 

NBAC3IR4NOLO 

CALL IOP(2HWB,a,LPUa, 2) 
CALL I0P(2HWB,4,NSUB, 62) 
NWS s NWS 4 64 
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GOTO 3 


2 NW aL0fMQfNE4« 


XF( IT.NE.NW ) 


IN SNW / 6a 


IFt MQ.EO, 9 ) NW| SNN • U 
NW2aNW«NWl*HQ<ia 

LPuasiRtLPua 

NBACaJR^KiOLO 

IF(NW«NWl*NW2»H0*a ,NE. 0 ) PRINT 
FORMAT( 6X *ERROR IN UNITR* ) 

CALL 

iFt Mo.GT. e ) cl!;[; 

tF( H0.6T. e ) CLL 

NWS I 

NC6aaIR 

NOLOalR 


lOP(2HWB»a,LPUfl, 4) 
I0P(2HWB,4,I0UH, NWl) 
lOP(2Hw0,a,B , MQ) 
I0P(2HWB,4,A , NN2) 
NWS ♦ NW 


RETURN 

END 

subroutine READ4 
CONMON 


COMMON 

common 

common 

common 

COMMON 


TRAC 

SCRl 

FRNT 

SU8T 

WRIT 


dimension 
EQUIVALENCE ( 


( ICO ) 

BC7M0J ,A(13P00) 

NWS,NWR,NBS,NUP 

IRf LPf II f I2f NR^HOf NW,S (2041 ) 

LD»MQ,ko,NQL,NQ^NQR,MP1 ,lQT(950) 

ISUB(^5) 

LPU4,NBAC, NWl, nW2,NC64, NOLO 
IDUH(956) 

lOUM ) 

COTO ( 1,2,4 ) 
CALL I0P(2HRB,4,LPU4, 2) 

CALL I0P(2HRB,4,NSUB, 62) 

NWR a NWR ♦ 64 


IF( MO.CT, 0) 

IF( MO.GT, 0) 

Nw a Nwl4Nw2fM044 


CALL lOP (?HRB,a,LPU4, 4) 
CALL I0P(2HRB,4,I0UM, NWl) 
call I0P(2MRB,4,8 , MQ) 

CALL I0P(2HRB,4,A , NW2) 

NWR a NWR 4 NW 

CALL I0P(2HSP,4,LPU4«NBAC) 

NBS a NBS ♦ NBAC 


CALL I0P(2HSP,4,LPU4«NC64) 
5 RETURN ♦ "'C64 

END 

subroutine store ( MTAPE,LPO8,K1,L1,0,A ) 
dimension B(1),A(1) j 

( «^OR / XR,LP,J1,J2,NR,MO,Nw 
common / TRAC / NWS, NWR, NBS, NUP 


IF( IT.NE.NW ) 


IR a NW / 64 

IT a IR * 64 

XR XR ♦ 1 
NW a XR *64 » 7 

NA a NW • nr 

LPOS a IR * LPOS 

LP a LPOS 
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CALL I0P(2Hwfl,NTAPe,lH,7) 

CALL IOP(2HW0,nTAPE,R(KI),NR) 
CALL I0P(2H«B,KTAPE,A(LI),NA) 
NWS 3 NWS ♦ NW ♦ 7 


RETURN 

ENP 

SU9R0UTINE SEMBl ( NTAPE, IGO,LPOS ) 

common f TRAC / NW8*NWR,NBS»NUP 

COMMON / SCRl / IR,I.P,Il,12,NR,MO,NW,S(20ai) 

GOTO ( U2»3 ) 
LP08 a IR ♦ LPOS 
LP 3 i.P08 

CALL I0P(2HWB,NTAPEf IR , 7) 

CALL 10P(2HwP,NTAPE,S , NW) 
NWS S NWS 4> NW ♦ 7 

CALL X0P(2HRB,nTAPE# IR » 7) 

CALL I0P(2HRB,NTAPE»S * NW) 
NWR a NWR ♦ NW ♦ 7 

call I0P(2HSP,NTAPE,LP-LP0S) 
N6S a NBS ♦ LPOS 


ICO 


GOTO U 


GOTO a 


U RETURN 
ENP 

subroutine SEMB2 ( NTAPE, IGO,LPOS ) 

common / TRAC / NWS»NWR,NB8,NUP 

COMMON / 8CR2 / IR,LP,Il,I2,NR,MO,NW,S(2a4in 

GOTO ( 1,2,3 ) XGO 


1 


2 


3 


CALL I0Pf2MWB,NTAPE,IR , 7) 
CALL lOPf2HWB,NTAPE,8 , Nw) 
NWS a NWS ♦ NW ♦ 7 

CALL I0P(2HRB,NTAPE,IR , 7) 

CALL IOP(2HR8,ntAPE,S , NW) 

NWR s NWR 4 NW ♦ 7 

CALL I0P(2HSP,NTAPE,LP*LP0S) 
NBS a NBS ♦ LPOS 


GOTO R 


GOTO a 


n RETURN 
END 

SUBROUTINE UNITED f NT APE , IGO, MD, B) 

DIMENSION BCD 

common / TRAC / NW8»NwR,nBS,NUP 

goto (1,2) ICO 

1 CALL XOP (2HMB,NTAPE,M0, 1) 

CALL top (2HWR,NTAPE#B ,M0) 

NWS a NWS A MO ♦ 1 

. GOTO 3 

2 call IOP (2HRB,NTAPE,ms# 1) 

CALL IOP (2HRB,NTAPE,B ,M8) 

NWR a NWR 4 MS A 1 


3 RETURN 
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ENO 

SUBROUTINE WIND ( NTAPE ) 

CALL lOP f SHREW, NTAPE) 

RETURN 

ENO 

subroutine PRNTED ( J 6 O,NUHB,M 0 I ,EO ) 
COMMON / PBOB / MAXF,MAXB,MAXW, 1 PRNT 
dimension EO(l) 

IF (IPRNT,EQ, 0 ) GOTO 1 

IF ( JGO.EO.l ) PRINT 100 , NUMB 
IF ( JG 0 ,E 0,2 ) PRINT 200 , NUMB 

PRINT 300, (EO (I), 1st , mod 

NUMB a numb » 1 


format ( 6X 

format ( 6X 

FORMAT ( IBX 

RETURN 

end 

subroutine PRNT0 

common 

common 

COMMON 

common / FRNT / 
common / PRNT f 
NUMl 8 NUMl ♦ 1 


common / frnt / lo,mo,ko,nq 
common / PRNT f IP0,IPI,IP2 
NUMl 8 NUMl ♦ 1 
IF ( NUMl.LE.NUM2 ) goto 99 
IF (IP0 ,EQ, 0) goto 99 
PRINT 10,LZ,NODEI,MD1 
PRINT 20,M1,k1,N1 
print 30, mo, ko, no 
format ( *1* ) 

format ( // 6X *L2 a* 19,6 

format ( 6X *Ml 8* ju,b 

FORMAT ( 6X *MO a* J9.6! 


*SUBSTR * 
♦ELEMENT* 
I 2 F 5.1 


( LZ ) 

NDOF (2500), 18 (2500), IC (2500), LQ( 150), LOC (150) 
NOOEl, 101(150), Ml, Kl,Nl,NOl,MOCl (150) 

N00E2, IO2(150),M2»K2,N2,NQ2,MOC2(150) 
LO,mo,KO,NQL,NQ,NOR,M01 ,LQT(950) 

IP 0 » IPI , IM 2 , IPA,IPB,NUM 1 ,NUM 2 , IPS 


10 format ( // 6X < 

20 FORMAT ( 6X i 

30 FORMAT ( 6X t 

99 RETURN 
END 

SUBROUTINE PRNTl 
common 


19 , 6 X *NO 0 El 8 * 19 , 6 X *M 0 l 
19 , 6 X *K 1 a* 19 , 6 X *N 1 


COMMON 

COMMON 

common 

common 


common / FRNT / LD,MQ,kO,NQ 
common / PRnT / IP0>IP1,XP2 
IF ( NUM1,LE,NUM2 ) GOTO 99 
IF (IPl ,E0, 0) GOTO 99 
PRINT 10,(101 (I),Ib1,N00E1) 
PRINT 20, (LO (I),Ib1,nO 0E1) 
PRINT 90,(MOCl(I),l8l,Nl ) 
PRINT 50,(MOC2(I),Ial,N2 ) 

PPINT 60, (LOC (I),I81,N1 ) 

format (// 6X *10 * 20I<| ) 

format ( 6X *LQ ♦ 2019 ) 

format ( 6X *M0C1* 2019 ) 

FORMAT ( 6X *M0C2* 2019 ) 

format ( 6X *LOC • 2019 ) 


19 , 6 X *KQ a* X 9 , 6 X *NQ 


NDOF (2500 ), IB ( 2500 ),IC( 2500 ) ,L0( 150), LOC (150) 
NODEl,lOl(150).Ml,Kl,Nl,NOl,ioCt(150) 

NOOE2, IQ2( 150) ,M2,K2,N2,NQ2,MOC2(150) 
LD,mq,kO,NOL,NQ,nGR,MO1,LOT(950) 

IP0» IPl, IP2, IPA, 1PB,NUM1,NUM2, IPS 


<w«nal paqc is 

Of fOOR QUALITY 
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\ 

10 

11 

99 


1 


a 

10 

20 

99 


100 

99 


SUBROUTINE PRNT2 ( LZ 5 

COMMON / FRNT / LD,MQ,KO,NQL»Nfi,NOR,HDl,LOTf950) 
common / PRNT / IP0,IPI,IP2,IPA,IPB,NUM1,NUM2,IPS 
IF ( NUMI,LE,NUM2 ) GOTO 99 
IF CIP2 ,EQ, 0) GOTO 99 
PRINT 1,LZ 

PRINT 10, (LOT (I), 1 = 1, MOD 
II = MOl ♦ I 
I? a MDl ♦ NQ 

PRINT ll,aOT(I),I = Il,12 ) 

FORMAT ( 6X *LZ * I« ) 

FORMAT ( 6X *LQT * 20I« ) 

FORMAT ( 6X *LOCT* 20ia ) 

RETURN 

END 

SUBROUTINE PRNTA (A,B,NQ,NT) 

COMMON / NOOO / NOtaSB) 

common / PRNT / IP0,IPl,IPH,IPA,IP0,NUMl,NUM2,IPS 
dimension A(i),en),o(i0B) 

IF (IPA ,EQ, 0) GOTO 99 
PRINT 20 
DO I 1=1, NT 

Qcn=0 

PRINT 10, (A(I) ,Isl,Nn,8(l) 

IF ( NQ.EO.l ) GOTO 99 

DO 2 K=2,NQ 

I1=K-J 

I2=NO(K) 

U=I2*NT-K 

PRINT 10, (Q(I) ,1 = 1,11), (ACn,I = I2,13),BCK) 

FORMAT ( 6X *A * 17F6,0 ) 

format ( / ) 

return 

END 

SUBROUTINE PRNTB ( B,NO ) 
dimension BCD 

common / PRNT / IP0,IP1,IP2,1PA,IP0,NUM1,NUH2,IPS 
IF ( IP8,EO,0 ) GOTO 99 
PRINT 100, (0(1), 1=1, NO) 

FORMAT ( / 6X *B * 10F6.1 ) 
return 


END 

SUBROUTINE PRNTS (NGO,LZ,N0, JPOS, IOEC,LP) 

COMMON f PRNT / IP0,IPI,IP2,IPA,IPB,NUM1,NUM2,IPS 
COMMON / FRNT / l.O,MQ,KQ,NQL ,N(l,NQR,MDl ,LQT (950) 
COMMON / STOR / IR,LN,J1,J2,NR,MD,NW 
common / 111! / R1,R2,RM,WM,H1,H2,NE,NC 
INTEGER R1 ,R2,RM,wM,HI,H2 


IF ( IPS,EQ,0 ) GOTO 6 

GOTO ( t,2,3,a,S ) N60 

1 IF ( N0,EO,1 ) PRINT 10,LZ,MQ,NQ,NE 

PRINT 15,NB,R1,R2,NC 

GOTO 6 

2 PRINT 30, LZ, JPOS 

GOTO 6 

3 IDEC X lOEC ♦ 1 

IF ( IDFC.EO.I ) PRINT 100, LP 

GOTO 6 

a PRINT 200, LP 

GOTO 6 

5 II = J1 ♦ MO 
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12 a J2 ♦ MQ 

print 11,12, 


If’ 

15 

30 

100 

200 

300 


fSrmIJ { ' IfJ ;J|b* {H'’; *"?* {n*’'' ‘"o* i«.jx .ne. u 

*“ *” ' ‘if-os* “ ‘"'=* 

poumJ! I 1!* *®EGINN1NG position IS* 16 ) 

format J III POSITION IS* 16 ) 

T ( 12 X *Il*I«»3X*l2*I<l,3X*IR*I«,3x#i.p*ia ) 


6 RETURN 
END 


APPENDIX C 
FORTRAN LISTINGS OF 
SKYLINE EQUATION SOLVERS 
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appendix c . 1 

SKYSOL SKYLINE EQUATION SOLVER 




subroutine SKYFaC (A,N0EG,nEND,N,UO,V,SINGA0, 

. PDCHEK,ICONO,lAFlLE,ACOND,OtTCF,IOETEX,NEGElC,IFAlL) 

C 

C TYPE AND DIMENSION STATEMENTS 

C 

REAL A(n,V(l) 

real aij 2 ,d,oetcf,ootpro,epsmac 

INTEGER LOCI) 

LOGICAL POCMEKrSINGAB 
equivalence (A1J2,D,DMAX) 

DATA EPSMAC/7.11E-15/ 

C 

c initialization 

C 

1FAIL=0 

IDETEXS0 

NEGEIGeP 

C 

C SAVE original matrix if IAFILE NE 0 and NBEG50 
C 

NNAsIABS(LO(N>n ) 

IF(IAFILE,EQ,0) go to 200 
IF(N0EG,NE,0) go to 200 
rewind IAFILE 

WRITE(IAFILE) (A(J),Jsl,NWA) 

C 

C COMPUTE SQUARED LENGTHS OF UNCONSTRAINED ROWS N0EG+I THRU N 
C 

200 N0EGP1=NBEG+1 

DO 1000 I=N0EGPt,N 

ii=Loci*n 

IF(II) 1000,1000^^00 
000 V(I)sA(II )**2 

MsII-1 

K3MAX0(NREGP1 , lABSaOd) )«M4-1) 
l=min0(neno, n-1 

IF(K-L) 500,500,1000 
500 DO 600 JsK,L 

IF(LD(JtD) 600,600,600 
600 AIJ2sA(M+J)**2 
V(I)=V(I)+AIJ2 
V(J)3V(J)4.AIJ2 
800 CONTINUE 
1000 CONTINUE 
C 

C FACTORIZATION SECTION 

C 

DO «000 JSNBEGP1,NEND 
C 

C COMPUTE KU SUPERDIAGONAL ENTRIES OF J-TH COLUMN OF tUl 

C IF UNCONSTRAINED 

C 

JJ3LD(J*1 ) 

IF(JJ) 41000, <i000, 1200 
1200 DsA(JJ) 

JMJalA0S(LD(J)) 

JKSJJ-JMJ 
KUsJK*! 

1F(KU,EO*0) go to 2200 
DO 2000 Ks1,kU 
laJ.jKfK 




✓ 


V(K)sB,Jl 

II=LD(Un 

IFCII) 18J10 

1000 M=HIN0CIl-iA0S(Lt)(i)),K)i<l 

IJsJMJ+K 

V(K)sAnj)-DOTPRO(ACII-M) ,V(K.M)-M) 
A(U)sV(K)*A(I 1) ' 

2000 CONTINUE 
C 

C COMPUTE DIAGONAL ELEMENT 0(J) 


0=0-D0TPR0(A(JMJfl),V,KU) 

w 

c singularity test 

c 

2200 T0LR0Ws8,*EPSMaC*SQRT(V(J)) 

IF(A0SCD).GT,TOLROw) go to 2500 
IF(SINGAB) go to 6000 
DsTOLRQW 

2500 A(JJ)si,0/o 

C 

c UPDATE determinant IF C.^TCF IS NONZERO 

3000 IF(DETCF,EO,0,0) go TO 3500 
DETCFsDETCF^O 

3200 IF(ABS(OETCF),LT,1,0) GO TO 3ft00 
OETCF=DETCF*,0625 
lOETEXslOETEX+a 
GO TO 3200 

3a00 IF ( ABS (DETCF ) ,GT , ,0625) GO TO 3500 
DETCF=DETCF+16, 
lOETEXslOETEX-a 
60 TO 3«00 

C 

C POSITIVE OEFINITNESS CHECK (IF PDCHEKs , TRUE , ) 

w 

350^ IF(O,GT,0,0) Go TO aMP0 

NEGEIGsNEGEIG-M 
IF(PDCHEK) go to 6500 

«000 CONTINUE 

C 

c SAVE factorization if IAFILE NE 0 AND NsNEND 


IF(1AFIlE,EQ,0) go to 5000 
IF(NEND,N£,N) GO TO 5000 
WRITE(IAFILE) (A(J) , j=i,nwa) 
IF(ICONO,EO,0) go TO 5000 

MATRIX CONDITION ESTIMATION 
KS0 

DHAXS0,0 

DO a200 Isl,ICONO 

SET UP random RHS VECTOR IN V 


original 

^ POOR 


IS 

quality 



CALL SkyrDM(V,NWA,0,0,k) 


SOLVE FOR V AND DO ONE 
OBTAIN 2»N0 RM relative 


CYCLE OF ITERATIVE 
SOLUTION ERROR IN 


refinement to 

DELTA 


CALL SKYSOL 


(A,N,LO»OOTPRO»^0^^»l,V,V(Ntl ),1 , lAFlLEr 






t 




I 


f 


( 


c 

c 

c 


, V(2*N>n,0FLTA) 

OMA?<sAMAXI (OELTA,OMAX) 

4200 continue 

NOW estimate C(A) from largest delta and machine precision 


ACONOsOHAX/( (1 ,0+DMAX)*EPSMAC) 

5000 RETURN 
C 

C ERROR EXITS 

C SINGULAR MATRIX(SlNGA8s,TRUE, ) 

6000 IFAILsJ 

OETCFs0,0 
60 TO 5000 

C INDEFINITE matrix (POCHEKs.TRUE, ) 

6500 IFAILs.J 

GO TO 5000 
END 

SUBROUTINE SKV80L(A,N,LD,OOTPRD,10P,IBX,B,X, 
• IREF,IAFILE,V, DELTA) 

C 

C TUYPE and DIMENSION STATEMENTS 

C 

INTEGER LDH) 

real A(n,B(l),BI,8XFAC»DELTA 
REAL RN0RM,V(1 ) ,X(1 ) ,XI,XNCRM 
EQUIVALENCE (BI,XI,XN0RM) 

C 

C INITIALIZATION 
C 

ASSIGN 3100 TO NEXT 

KREFS0 

BXFACS0.0 

IF(lBX*EQ,0) GO TO 200 
BXFAC=1,0 
DO 150 Isl,N 
150 X(IJ=8Cn 

200 IF(IOP,GT,0) GO TO 1000 
IF(IBX.EQ.0) go TO 1100 
C 

C RMS MODIFICATION 
C 

DO 1000 Isl^N 

ii«Lon«i) 

IFCIl) 300» 1000, 1000 
300 BIsbCD 

IF(BI.EQ,0,0) GO TO 1000 
Ils-ll 

KsI-II+IABS(L0(I))+1 
DO 900 JsKfN 
JJeLD(J^l) 

IF(JJ) 900,900,400 
400 MsJ-l 

IF(M) 900,600,600 
500 X(J)aX(J)«»A(IJ+M)*Br 
GO TO 900 
600 IJsJJ-M 

IF(IJ-IABS(LD(J))) 900,900,000 
000 X(J)aXCJ)-A(TJ)*flI 
900 continue 

1000 continue 






c 

c 


forward substitution pass 


1100 


1300 

1300 


1500 


C 

C 

C 


DO 1500 IslfN 

ii=LD(r+n 

IF(II) 1300,1300,1300 

X(I)S0,0 
60 TO 1500 
IMIslABSfLD(D) 

McII-imI-1 

J”J=J«JJ’D0TPRD(AnMl4l),)((i.H),H) 

IF(1OP,NE,0) go to 5000 

SCALING PASS 


1000 00 3000 Isl,N 

IlalABSaOd + D) 

3000 X(I)sA(II)*x(l) 

C 

C BACKSUBSTITUTION PASS 
C 


IsN 


DO 3000 Ksl,N 
ii=LO(i+n 

IFCIl) 3300, 3300, 2U00 
3300 X(I)s0XFAC*B( 1) 

GO TO 3600 


C 

C 

C 


3«00 MsII-t48S(LD(I))-l 
IF(M,eo,0) go to 3000 

00 3500 Jsl,M 


X(I-J)=X(I-J)-A(II-J)nx(n 

3500 continue 
3800 IsUl 
3000 CONTINUE 

IF(1BX*IAFILE,EO,0) go to U000 

GO TO NEXT, (5100,3500) 


iterative refinement section 


3100 


C 

C 

c 

c 


KREF=KREF+1 

IF(KREF-IreF) 3300, 3300, <1000 

calculate RESIDUAL VECTOR tR)0s[Bl-(A] 
TRl RETURNS IN X, AND OLD IXJ IN V 


IX) USING 8KVMUL* 


C 

C 

C 

C 


3300 rewind IAFILE 

NWAsIaB5(L0(N+D) 

READ(IAFIlE) (A(J),Js1,NWA) 

CALL SKVHUL (A,N,L0,X,V,-1,B,V) 


SOLVE Fos CORRECTION IPX! (WHICH AFPEXRS IN X). CORRECT 

OLD SOLUTION AND EVALUATE 2-NORM RELATIVE ERROR DELTA 


REAO(IAFILE) (A(J) , Js 1 ,NWA) 
BFAC«0,0 

assign 3500 TO NEXT 
CO TO 1100 
3500 XNORM80.0 
RNORH>0,0 
DO 3600 lsl,N 
RNORM»RNDRMfX(I)*X(I) 


c 

c 

c 


' x(i)=vrn+x(n 

XNORMsXNORMfX(I)*X(I) 

3R0« CONTINUE 
OELTAse.pl 

if(xnorm,ne,0,0) del t Assort (rnorm/xnorm) 

IF(KREF-a) 3100. 3l00»«f500 


CONSTRAINED RHS RECOVERY 


«000 IF(IBX.LE,0) GO TO 5000 
DO 0000 Isl.N 
II=LD(U1) 

IF(in 0200.0200.0000 
0200 IHIslARS(LOtl) ) 
Ms-IUimi-1 


B(I)sDOTPRO(AnMI + l).x(I-M),H) 

DO 0600 Jsl.N 
IJslA0S(LD(J*l))+j-l 
IF(IJ-A0SfLD(J) ) ) 0600.0600.0000 
0000 P(I)sH(I)fA(IJ)*x(J) 

0600 CONTINUE 
OP00 CONTINUE 
5000 RETURN 


END 

subroutine SKYMul (A.N.LD.X.AX.IOP.BfR) 
double PRECISION AXm.AIJ.AXl 
INTEGER LOm 
real Am. 6 (i).x(n.R(i) 

DO 2000 Isl.N 


ir = iARS(LD(i + n) 
AXmsDBLE(A(in )*OBLE(X(I) ) 
M=II-IABS(LD(n )-l 
IF(M,EO,0) go TO 2000 
DO 1500 KsI.m 


Js I 


aIJ=A(II-k) 

AX(I)=Axn)fAIJ*OBLEfX(J)) 

AX(J)sAX(J)+AIJ*DBLE(X(I)) 

1500 continue 
2000 CONTINUE 

IF(IOP) 2500»5000. 3500 
2500 DO 2800 Isl.N 
AXIsaX(I) 

Rn)sx(i) 

X(naOBLE(BCm-AXI 
IF(LOtI+1),LE.0) X(I)s0, 

2800 CONTINUE 

GO TO 5000 
3500 DO 3000 Isl.N 

R(I)sDBLE(Bfin-AX(I) 

IF(LD(I+1 ),LE,0) 0(1)30, 

3800 CONTINUE 
5000 RETURN 
END 

subroutine SKYROM (V.N.VMEAN.IROM) 

real V(N) 

C*VmeaN- 0,500 
IFdROM.NE.B) GO TO 300 
IXS26903 
IROMsl 

300 00 500 Isl.N 

iXsMOn (27 181 *lXi 13089.65536) 
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XsIX 

V(I)a,l5250789E-«*XfC 
500 CONTINUE 
RETURN 
END 


APPENDIX C . 2 

GASP COLUMN EQUATION SOLVER 
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PROGRAM MAlNl ( INPU1#0UTPUT»TAPE1»TAPE2»TAPE3»TAPE<I 

#TAPES,TAPE6,TAPE7, TAPES ) 

COMMON NOOF(ia«n 

COMMON / A / i0(a),N00ES,jQ(a),s(e,6) 

COMMON / 8UPORT / I60(«) 


00 to I=i#B 

iO CALL ZOQIN ( bHREWINOft ) 

READ 10O»NOFRE, nodes, NOIVZ rNOIVYfNOlVX.lOiNOSP 
noeu*noivy*noivx 
NOPT a<NOIVY + n*(NOlVXft ) 

PRINT 20O»NOPT,NOFRE»NOEL»NODES»NOIVZ,NOIVY,NOIVX,IO,NO3P 

CALL I0BIN(5HWRITE»6,N0PT,l) 

CALL I0BIN(SHMRITE,6«NQEL, t) 
call X0BIN(5HwRITE*6,N0SP»1) 

DO I Isi,NOPT 

1 NU0F(I)3NDFRE 

CALL IO»IN(5HwRITE,6,NOOF,noPT) 

READ 100, JO 
print 300, jo 
NP laNOIVY^l 
DO 3 lXal,NOIVX 

JQ(na(ix-n*NPi+i 
Jo(«)=jo(n + i 
jQ(2)aJW(n>NPJ 
JO(3)ajO(a)*NPl 
DO 3 IY5l,wplVY 
00 2 Isl,a 

2 lQ(I)sJQ(I)^IYi*l 

call IO0lN(5HwpiTEf5. IQ,«) 

3 continue 

M a N00ES*N0FRE 

00 5 Ixt,M 
DO 5 Jal,M 

5 S(I,J)8l ,0E»20 

00 6 lat ,M 

S 8(I,I)at,0 

00 7 L2sl,N0EL 

7 call lOaiN ( 5HNRlTE,6,S,6a ) 

DO 6 tai,NOSP 

read too, IBO 

8 call I08IN ( 5HWRlTE,6,IB0,a ) 

CALL SECOND ( T ) 

call JOBInFo(8,J1) 

CALL GASPl ( 10 ) 

CALL JO0INFOC8. J2) 

CPaFLOAT(J2.Jl)/l000, WffPfWMkiY 
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CALL TINEX ( 0,0, T ) 

PRINT a00,T,CP 

<♦00 FORMAK 6X *TXKE FOP GASP • F6.2/7X * CP FOR GASP * F6,2) 


IBM FORMAT 

( 

1013 ) 




200 format 

( 

* 1 * 

5X 

*nopt * 

15/ 


• 



6X 

*NOFRE* 

15/ 


t 



6X 

♦NOEL ♦ 

15/ 


0 



6X 

•NODES* 

15/ 


0 



6X 

*N0IVZ* 

15/ 


0 



6X 

♦NOIVY* 

15/ 


0 



6X 

♦NDIVX* 

15/ 


0 



6X 

*10 ♦ 

15/ 


0 



6X 

*NOSP ♦ 

15 

) 

300 FORMAT 

( 

/// 

6X 

♦JQ ♦ 

6X5 

) 


END 

subroutine get ( 01»02»G3»U4,LQ»LO»LI*NUM ) 
DIHENSION LQ(5000)*LO(5O00) 
common / G / IGET(26) 

INTEGER U1,q2,q3,Q<i 

CALL XOOIN ( 4HREA0#4, IGETr26 ) 

OlsXGETCl) 

Q2alCET(2) 

03aIGET(l) 

OasIGETCa) 

LI B NUH^t 
NUMsNUH * 03 

Xl«4 

X2»I1^03 

00 1 L=LI»NUh 

IlBllfl 

l2Bl2fl 

LQ(L)«IGET(I1) 
t LD(L)sigET(I2) 

return 

END 

subroutine TIMEX ( U»NUM,Tt ) 

IF ( NUM.EO^B ) NUM s 1 
call second ( T2 ) 

Tl B T2 - T1 

IF ( lI.EQ.R ) goto 10 

FAC B num 

FAC b ( Tl / FAC ) * 1,0E6 
PRINT lPi3,II,NUM,Tl,FAC 
100 FORMAT ( 6X *H * IS • 


t 


2X *NUM • X5 , 

2K *TlhE* F9,2, 

2X *FAC * FIS, 2 ) 


10 RETURN 
ENO 

subroutine GA8P1 ? IC ) 

COHHON XB(iaai),JBnoai),NOOF(iafln.l8PT(ia<»nrR(«323) 

COMMON / A ✓ 20(R)»NLOE8f J0(<1),8(8,8) 

If ( IO,GT,0 ) PRINT 200 

CALL MIND ( 5 ) 

CALL mind ( 6 ) 


CALL lOBXN ( 
CALL lOBIN ( 
call lOBlN ( 
XF ( XO,GT,0 


«HPEA0,6,N0PT,1 J 
RHREAO,6,NOEL»1 ) 
RHREAD#6,N0SP, I ) 

) PRINT 210,NOPT,NOEL*NOSP 


CALL lOBIN ( RHREAO»6,NOOF»NOPr ) 

IF ( IO,6T,0 ) PRINT 200 

IF ( IO.GT .0 ) PRINT 220,(NOOF(n,Isl,NOPT) 

NOEQ s 0 


00 2 I«1,N0PT 

IB ( I ) « 0 
IBPT f I ) s 0 


NOEQ » noEO ♦ NOOT ( X ) 


IF ( IO,GT,0 ) PRINT 230, NOEQ 
IF ( IO,Gr ,0 ) PRINT 200 

00 « LZ»1,N0EL 


CALL lOBIN ( aHPEAO,5,IO,<» ) 

IF ( IO,GT,0 ) print 2R0,LZ»NOOE8»(IO(I),I«j,noOES) 
00 3 JJ«t#N00ES 


J ■ in ( jj ) 

IS f J ) » IB ( J ) ♦ j 

00 3 xisl, NODES 

I ■ 10 ( II ) 

IF ( I,LT,J ) GOTO 3 
L a ( I»J ) ♦ i 

IF ( L.CT.IBPTrj) ) I3PT ( J ) . L 

3 continue 

R CONTINUE 

J1 ■ t 

00 6 I«1,N0PT 


J2 « J1 ♦ NOOF ( ( ) •! 
X * IB ( I ) 

DO 5 JaJj,J 2 
5 R ( J ) s X 

*» Jl a J2 ♦ 1 


CALL lOBlN; ( 5mwr1TE>2#H»N0EQ ) 

PRIWT 2(?l0 

PRINT 250#(IB (X)f Isl^NOPT) 
PRINT 200 

PRIM 260,(IBPT(I),Ial,NOPT) 
print 200 

call wind ( 5 ) 

CALL WIND ( 2 ) 


IP ( lO.CT.O ) 
IP C XO,GT,0 ) 
IP ( IO,GT,0 ) 
IP ( IO.GT.0 ) 
IP ( lO.GTfB ) 


CALL OECON ( NOPT,IT,ro ) 

IP ( IO,CT,0 ) PRINT 270,11 


wc;t,ur>u i I 1 ] 

call JO02NFO(8,Jl) 

CALL SEARCH ( M0EL#N0PT,U» 
CALL JO8INFO(0,J2) 
SCPHsFLOAT(J 2 «jn/l 000 , 

CALL TIMEX ( 0,0, 


Tl ) 


CALL SECOND ( T2 ) 
call J 0 BlNF 0 C 8 ,jn 
CALL FORHQO ( NOEL#NOPT, 
call JOBINFO( 8 , J 2 ) 

PCPOaFLOAT(j2.ji)/j(5p0^ 

call TIMEX ( 0,0, T2 


10 ) 




call second ( T3 ) 
call J0BInP 0(8, Jl) 

CALL FORHdS r NOSP, IT, 10 

Call jobinFv*8,j 2) 
PCPSaFLOAT(j 2 «jn/l 000 , 

CALL Timex t 0 , 0 , T3 ) 


CALL second ( T« ) 

CALL JOBInFo( 8 , Jl) 

CALL Reduce t nopt#noeq,it,io 

CALL JOBINFO( 8 ,J 2 ) 
RCPEsFLOAT(J 2 «Ji)/UJt»o, 

call TIMEX ( 0,0, Tfl ) 


^ CALL SECOND ( T5 ) 
CALL JOBINFO( 0 ,Ji) 

call Solve ( nopt,noeq, 
CALL J0BINF0(8, J2) 
8 CPE«PLOAT(j 2 -jn/ 1000 , 
CALL TIMEX ( 0,0, T5 


10 ) 


) 


print 300,T1,SCPH 


original 

^ POOR 


■^GE »S 
QUALny 


PRINT a«0,T2,FCPO 
PRINT 500fT5#FCPS 
PRINT fc00,TO,RCPE 
PRINT 000,T5»SCPE 


I 300 Format ( 6X *time for search* f6,2 / 

ll , • 6X * CP for search* F6,2) 

f 400 format ( 6X *TlMe FOR FORMQO* F6,2 / 

!! • 6X * CP FOR FORHaO* F6,2 ) 

t 500 FORMAT ( 6X *TIM£ FOR FORM03* F6,2 / 

' • 6X * CP FOR FOPMOS* F6,2 ) 

600 format ( 6X *TIME FOR REDUCE* F6,2 / 

• 6X * CP FOR REDUCE* F6,2 ) 

600 format ( 6X *T1ME FOR SOLVE * F6,2 / 

• 6X * CP FOR SOLVE * F6,2 ) 


200 FORMAT 
210 FORMAT 


220 FORMAT 
230 format 
240 FORMAT 
250 FORMAT 

260 format 

270 FORMAT 


I *1* 10X **,.,,,GASPl.,.,,,* ) 

( 6X *N0PT* 15/ 

6X *NOEL* 15/ 

6X *NOSP* 15 ) 

( 6X *NDOF* 3013 ) 

( 6X *NOEQ* 15 ) 

( bX *tZ* 15, 5X *NOOES* I5,5X *10* 2013, 5X ) 
( 6X *IB * 3015 ) 

( 6X *I8PT* 3013 ) 

( 6X *IT* IS ) 


RETURN 

END 

SU8R0UTINE DECOM ( NOPT, IT, TO ) 

COMMON 19(1441 ),J8 (1441) ,NOOF (1441) ,I8PT( 1441) 

COMMON / R / Il,I2,Kl,K2,M0 
dimension 1C(5) 
equivalence ( 11, IC ) 

IP s 1 
IT s I 

11 s 1 

12 * IBPT ( 1 ) 

1 L a IP t IQPT ( IP ) - 1 

IF ( L.LTtNOPT ) CALL UPOOF ( L, IP,K1 ,K2 ,mo ) 

CALL lOBIN ( 5HWRITE,3,IC,5 ) 


IF ( IO,GT,0 ) PRINT 10,11,I2,K1,K2,MO 
IF ( L,GE,N0PT ) GOTO 4 

ICHEK c 70PT ( IP ) 


2 lU B IBPT(IP*1)-ICHEK 


IF 

( 

IU.GE.0 

ICHEK 

B 

ICHEK - 

IP 

B 

IP ♦ 1 

12 

GOTO 2 

~ 

12 - 1 

3 11 

8 

12 

12 

B 

12 * lU 


) 

1 


GOTO 3 


m 


GOTO 1 

IT 

3 IT 

♦ 1 


10 FORMAT 

( 

6X 

★ II* 

I6» 

• 


2X 

*12* 

16, 

• 


2X 

*K1* 

16, 

• 


2X 

*K2* 

16, 

«*RETURN 

end 


2X 

*M0* 

16 ) 

subroutine 

UPOOF 

( L» 

IP,Kl,K2,M0 ) 


COMMON 18(ia41),JB(l<i<il)»ND0F(ia(ll),lBPT(14ai) 

ICHEK s I9PT ( IP ) 

ID = IP 

Kl s NDOF ( ID ) ♦ I 
K2 = 0 
MO = I 
00 I J=IO,L 

J K2 s K2 ♦ NOOF ( J ) 

2 lU = I0PT C ID+l ) - ICMEK 

IF ( IU.GE.0 ) GOTO 3 

ICHEK s ICHEK - 1 

ID 3 ID ♦ 1 

Kl 3 Kl f NOOF ( 10 ) 

MO s MO > 1 

GOTO 2 

3 RETURN 

end 

SUBROUTINE search (NOEL »N0PT, I3AN0P, 10 ) 

COMMON 1 0( uan ,jB( I ^ai) , NOOF (i<iun , lq(i I » 5i ) , ld(i I ) 
COMMON /A/ lQ(4),N00ES>JQ(ti)»S(ef8) 

1 J31,N0PT 
J0(J) 3 I8(J) 

IF ( lO.GT.l ) PRINT 200 
FORMAT (IHl 5X *SEARCM*) 

13 IsUIBANOP 

LO ( I ) 3 0 

13 Jsl,5l 
LQ(I*J)s0 
IK 5 I 

n LZsliNOEL 

CALL lOBlN ( 4HREADr5» IQ,4 ) 

3 II3UN00ES 
I 3IQ(II) 

10(1) 3 18(1) - 1 

3 JJ31, nodes 

J 3 IQ(JJ) >. IK ♦ I 
IF ( IQ(JJ),GT,I) goto 5 

2 LL»lrIBAN0P 

IF ( LQ(LLrJ).EQ,0 ) LQ(LL#J)3l 
IF ( LQ(LL»J).EQ«I ) GOTO 3 
2 CONTINUE 


00 

I 

200 

DO 

DO 

13 

DO 

DO 

DO 

DO 


3 CONTINUE 

PRINT*3 i ^ CALL ORDER ( IK » I BANOP, 10, NOPT 

379 F0RMAT(5X,*lZ*,T6) 

U CONTINUE 

CALL WIND ( (I ) 

CALL '^INO f S ) 

RETURN 

end 

SUBROUTINE ORDER ( IK, IBANOP, I0,N0PT ) 

COMMON IB(l««n,ja(iaai),N00Fn9«l),L0(ll,5n,L0(ll) 

COMMON / G / IGETCa6) 

9 DO 5 lst,IBANDP 

ZF ( L0(!,l),EO,e ) GOTO 6 

5 IT s 1 

6 DO 8 Ial,iT 

BIN s LO(I,l) 

IOC = I 
00 7 JsI,IT 

IF ( LO(J, 1) ,6T,MIN) goto 7 
MIN s LQ(J,U 
LOG n J 

7 CONTINUE 

LO(LOC,n = LQ(I,n 
LQ( 1,1) s MIN 

8 CONTINUE 

Locn = 1 

IF ( IT.EO.l ) GOTO 19 
DO 9 L=2,IT 

I = LQ(L-1,J) 

9 LD(L) a LO(L-l) ♦ NDOF ( I ) anOOF ( IK) 

19 la LQ(IT,1) 

M0= LD(IT) ♦ N00F(I)aN00F(IK) - I 

IF ( lO.GT.l ) PRINT IJ00,IT,MO 

IF ( 10,GT,1 ) PRINT 1200,(LQ(l,n,lsl,IT) 

IF ( lO.GT.l ) PRINT 1250, (LOU), Isl, IT) 

IGET(1)sJ0(ik) 

IGET(2)sN0OF(1K) 

IGET(3)alT 

IGET(a)sMo 

Ils« 

I2SI1+IT 

DO 90 Ial,IT 

IGET(UI1)=L0(I,1) 

90 ICET(ItI2)sL0(I ) 

CALL lOBlN ( 5 HwrITE,9,ICET,26 ) 

IF ( IK.EO.NOPT ) GOTO 12 
DO 10 j3l,5B 
DO x0 Ial,XBANDP 
10 LO(I,J) a LO(I,Jfl) 

IK a IK ♦ I 

IF ( IB(1K),EQ.0 ) GOTO a 

1300 format ( 5X AMQ* 15 / 

• 5X *MD* 15 ) 


1200 format (SX *LQ* 16X5) 

1250 format (5X *L0* 16X5) 

12 RETURN 
ENO 

SUBROUTINE FORMQD (NOELf NQPT, 10 ) 

COMMON 18 (335)#NOOF(335),HQ(33S),HD(335)»LQ(5000),LD(5000)» 
• 0(37000) 

COMMON /A/ XQ(a),N0DESf JQ(6)»S(6r8) 

INTEGER OUM 

IF ( IO,GT,0 ) PRINT 200 
200 FORMAT (IHl 5X ♦FORMOO^k) 

IK « 0 
NUM S 0 
OUM a 0 
JOUT a 1 
00 R LZalfNOEU 

CAI,L lOBXN ( UHREaO»5,IQ»R ) 

CALL lOaiN ( 4HREA0f6rS >64 ) 

J = IQ(1) - IK 

00 6 JJal, NODES 

K a IQ ( JJ ) • IK 
IF ( K.GT.J ) J s K 

6 CONTINUE 

IF (J.GE.JOUT) CALL EXPAND (IK, J, JCUT>NUM,DUH) 

IF ( IO,GT,0 ) PRINT 1149, LZ, NUM, DUM 

DO 7 JJal , nodes 

K a 10 ( JJ ) - IK 

7 JQ ( JJ ) a NDOF ( K ) 


00 1 11=1, NODES 

I 3 lO(II) •• IK 
IB(X) a XB(I) • 1 
DO 1 jjal , nodes 

J = IQ(JJ) - IK 
IF (J.GT,I) GOTO 1 
CALL FIND1(X,J,L0Q,L00) 
CALL ADOS (II,JJ,LOD) 

1 CONTINUE 


2 IF ( XB(l),GT.0 .OR. IK.EQ.NOPT ) GOTO 4 
LT a HD ( 1 ) 

CALL lOBlN ( 5HWRITE,7,D,LT ) 

CALL SHRINK ( IK, I, JOUT, NUM, OUM ) 


GOTO 2 


4 CONTINUE 

call wind ( 3 ) 
CALL WIND ( 4 ) 
CALL WIND ( 7 ) 


1AA 




1109 


1 


2 

3 


1 

2 

3 

90 


1 

2 


3 


FORMAT ( 5X *LZ * IS, 

• SX *NUM * IS# 

, 5X *DUM * IS ) 

RETURN 

END 

subroutine expand C iKf J, J0UT,NUM,0 Um ) 

COM -ION IB(335),NOOF(335),MG(335),MO(335),LO(5000),LD(S000), 
, 0(37000) 

Integer dum 
00 3 n=jout,j 

CALU get ( IB(N) ,NOOF(n) ,M0(N) ,H0(N) ,LQ,LD,L1 »NUM ) 

DO I LsLl,NUM 

LQ(L) s LO(U) - IK 

LI a DUM ♦ 1 

OUM a OUM ♦ M0(N) 

DO 2 L=LI,0UM 
D(L) 30,0 
CONTINUE 
JOUT 3 J ♦ I 

return 

END 

SUBROUTINE FINOU I, J, LOO, LOO ) 

COMMON IB(335),NOOF(335),MO(335),MO(335),LO(S000),LDC5000), 
. 0(37000) 

LOO = 0 

LOO a 0 

L 3 0 

IF (J,EQ,i) GOTO a 
M2 a J*l 
DO I Mai, M2 

L 3 L ♦ hQ(M) 

LOO 3 1,00 + MO(M) 

LI 3 L ♦ 1 
La s L ♦ MQ(J) 

DO 3 LL=L1,L2 

IF ( LO(LL),EO,I ) LOO s lL 
CONTINUE 

LOO a LO(LOO) ♦ LOO 

IF ( LOQ.EO.0 ) PRINT 90, I, J, LOO, LOO 

format ( 5X *ERROR IN FIND# / 5X,«15/5X,2I5) 

RETURN 

END 

subroutine ados ( II,JJ,LOO ) 

common I0(335),NOOF(335),MQ(335),MD(335),LQ(5000),LD(5000), 
, 0(37000) 

COMMON /A/ IO(a),NODES, JQ(a),S(8,8) 

Ml 3 1 

IF (lI.EQ.l) GOTO a 
12 » II - I 
DO 1 i>i,ia 

Ml a Ml ♦ JQ(I) LWiGiNHi. = 3E IS 

M2 8 Ml ♦ JO(II) - I POOP 

N1 a 1 

IF (JJ,EQ.l) GOTO a 
J2 8 JJ . 1 
00 3 J31,J2 
Ml 8 N1 t JQ(J) 
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« N2 » Nl ♦ JO(JJ) • I 
00 5 

00 5 NsNl,N2 

O(tOO) a O(LOO) ♦ S(M,M) 
5 LOO a LOO f 1 

RETURN 
end 


SUaROUTXNE shrink ( IK|M0, JOUT,NUM,OUM ) 
INTEGER OUH 


XT a MO ( 1 ) ♦ 1 

LT a MO ( 1 ) ♦ 1 

K a 0 

J1 * HO ♦ 1 

J2 a JOUT • 1 

00 2 Jajl,j 2 

K a K ♦ 1 
I 6 (K) 3 XB(J) 

NOOF(K) 3 NOOF(J) 

HO(K) a MQ(J) 

2 MO(K) 3 MO(J) 

K 3 0 

00 3 Ja IT,NUM 
K a 1 

LOCK) a LO(J) - HO 

3 LO(K) 3 LO(J) 

K a 0 

00 U JaLT,OUM 
K a K ♦ t 

^ 0(K) a D(J) 

JOUT a JOUT • MO 
NUH a NUM • IT ♦! 

DUM a oUm • LT ♦ I 
IK B IK ♦ MO 
RETURN 
END 

subroutine F0RM03 < NOSP,IT#IO) 

”JJJ2j,N®'>^(«5).HO(J3S),HB(S35).UO(Sl)00),LI){S0Ba) 

COMMON /R/ n,I2,Kl,K2,M0 
COMMON /SUPORT/ IB0(«) 
dimension IC(5) 
equivalence ( n,ic ) 

Integer oum.row^run 

IF ( I0«GT,0 ) print 2000 

IK a 0 
NUM 3 0 
DUM 8 0 
ROW a 0 

CALL^IOBIN ( THREAD, 6, ISO, fl ) 

DO I lRal,lT 

CALL IOBIn ( ttHREAD r3#IC,5 ) 


fLD(S000), 


>0(37000) , 




CALL lOBlN ( SHWRITEr I , IC»5 ) 


call Explod (Ik,num,oum) 

IF C IO,GT,« ) PRINT 30«0,IR,num,OUM 

CALL STORE (NOSP,RUN,row) 

IF ( IR.EO.IT ) goto i 

CALL COMPACT (NUM.OUM) 

1 CONTINUE 

call wind ( 1 ) 

2000 FORMAT ( *1 FoRMDS * ) 

J0(.a FORM.T ( SX .IR. lt,2x »NUM» 16, 2X .OUH. 16 , 


RETURN 

END 


subroutine EXploD ( IK#NUM,DUM ) 

COMMON /R/ H,I2,Ki,K2 
INTEGER OUM 
DO 3 Nail, 12 


MO 


» 


CALL get ( IB»NO0F(N),Mi3Cn),M0(N),L0,LD,L1,NUM ) 

DO 1 LaLl»NUM 
1 UO(L) s LO(U) - IK 

L2 = DUM ♦ 1 
OUM = OUM ♦ MO(N) 

LIMaOUM»L2^1 

call IOBIn ( ^HREAO07,OCL2),LIM ) 

3 CONTINUE 

IK s Ik > MO 

RETURN 

End 


1 


2 

3 

4 


SUBROUTINE STORE 
common NoOF(335) 
• DS(100) 


(NOSP,RUN,ROW) 
»MQ(335) ,M0(335) 


COMMON / SUPORT / IB0(4) 
COMMON /R/ I1,I2,K1,K2,M0 
INTEGER ROW, RUN 
00 9 1=11,12 

MBANO s 0 
00 1 j=l,l 


» IQ ( 5000 ),LD (5000), 0(37000), 


00 


00 


MBAND B MBANO ♦ NOOF(J) 
2 J=1,I 
J1 3 J 

CALL FIN02(I,J,lOO,LOO) 
IF ( LOO,CT,0 ) GOTO 3 
MBANO 3 MBANO • NOOF(J) 
NWOROS 3 M04NO*NDOF(n 
^ K=I, NWOROS 
08(K) s 0,0 
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Ml a 

1 


00 

6 

Ja J1 

,l 




CALL 

FIN02( I,J, 

LOO, LOO ) 



IF ( 

LOQ.EQ.0 ) 

GOTO 6 



M2 a 

Ml ♦ NWOROS 

• M0AND 



NJ a 

NOOF(J) • 1 


00 

5 

MaMl, 

M2*MBAN0 




N2 a 

H ♦ NJ 


DO 

5 

Nari,N2 




OS ( 

M ) a OS ( 

N ) ♦ 0 ( LOO ) 



LOO 

a LOO ♦ 

1 



Ml a 

Ml ♦ NOOF ( 

J ) 



ROW a 

ROW f 1 

IF 

( 

ROW.NE.lBDd) ) 1 

GOTO 8 


Ml a MBANO •• NDOF(I) 

M2 a NWOROS ■» NOOFCI) 

LOC a 1 

00 7 M a M1 »m2#MBAND 

IF ( IBDCLOC^’l) •GT.B ) DS(M+LOC) a I.E50 
7 LOC a LOC ♦ I 

run's'*ru 5'‘I'?‘’*^ ’ ‘ «HRUD,6,ieo,« ) 


* CALL 1O0IN ( 5HWRITE»1,ND0F(I),1 ) 

CALL I08IN ( 5HWR1TE# I , NWOROS .1 ) 

CALL lOBXN ( 5 HWRItE,1,0S .NKCROS ) 

9 CONTINUE 
return 
End 

SUBROUTINE FIN02( I,J,LOQ,LOO ) 

LOQ a 0 
LOO a 0 
L a 0 

IF (J,EQ,1) GOTO 2 
M2 a J,l 
DO 1 Md1,H2 

L a L T HQ(M) 

1 LOO a LOO ♦ MOfM) 

2 LI a L ♦ 1 

L2 a L t MQ(j) 

00 3 LLaLtrLZ 

, IF ( LOfLD.EO.I ) LOO a lL 

3 CONTINUE 

LOO a LO(LOQ) ♦ LOD 

return 

end 

subroutine compact ( NUM,OUM ) 

’common /R/ IWI2,KWK2,M0 
Integer dum 

it a 1 

LT a t 
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DO 1 Jsl,MO 

ITa IT 4 MO(J) 
1 UTs LT 4 HO(J) 

K a 0 

Jl a MO 4 1 
J2 a la 


DO 2 J=JI,J2 

K a K 4 I 
NDOF(K) b NDOF(J) 
HQ(K) s MOfJ) 

2 MO(K) a MO(J) 

K T u 

Du i ja 

K a K4 1 

LO(K) a tOCj) - MO 

3 LD(K) s LOCJ' 

Kao 

DO it JaLT#DUM 
K a K 4 1 

« D(K) a 0(J) 

»^UM a mum • IT 4l 
DUM B DUM • LT 4 I 
RETURN 

end 


SUBROUTINE REDUCE ( NOPT# NOEO* I T» 10 ) 

COMMON 0(33000), jD(30«),KO(Jib0),M8(3tJP),NUM 
common /R / I1,I2,K1,k2,MP ' 

DSn«0),NOOF,MOEL.IB.LCPJ(K.«l) 
equivalence ( II, ic ) 


LOC a 0 


MI a I 
LI a J 

DO it lRal,iT 


NCALaO 


CALL SECOND ( 01 ) 


CALL IOBIN ( RHREA0,1,IC»5 ) 


IF C IO.GT.0 ) PRINT 200,IR,I|,12 


Do 3 11 * 11,12 

IS a KD(MO) 

ILS a JO(MD) 

XLT a J0(M2) 

IE 8 •JD(lS4l) - IS 4 KO(IS+l) 


IF 


IB 4 NOOF ,NE, 3 ) GOTO 10 
0 ( ILS ) a 0 ( ILS ) / 0 ( ILS- 
D ( ILSTI) a SORT ( 0 ( ILS 4 1 ) 


1 


) 

0 ( ILS )**2 ) 


NUM a 1 
GOTO 2 


10 


NUH B 0 

IF ( XB.GT.NOOF ) GOTO 1 
NUM 8 0 

D(ILS) a SORT ( 0(IL8) ) 

IF (MOOF.EQ. 1) GOTO a 

I * ILS ♦ NOOF 
P(X) a D(I) / D(ILS) 

Od^n a SORT ( Odfl) - Od)**? ) 
XLS a ILS ♦ 2*Nr>0F 
L s L ♦ 2 

NUH 8 NUH ♦ 1 

IF (NOOF,EO,2) GOTO 2 

1 CALL SECOND (Tl) 

CALL CHOLa ( L,IS#XLS,lLT,Xa»IE ) 

C CALL TXMEX ( XI»Nuh,t1 ) 

2 NCAL a NCAL ♦ HUH 

LOC a LOC ♦ 1 
L ad 


DO 20 NRLn»L2 
L a L ♦ I 

20 D8 ( L ) a 0 ( N ) 

CALL lOBIN ( bHwKlTER,0.OS, 10J,LCFJ ( LOC ) ) 


3 CONTINUE 

IF ( IO,GT,0 ) PRINT 300, L2 

IF ( IR,LT,IT ) CALL UPDATE 
CALL TIMEX ( IR,NCAL,QI ) 

0 CONTINUE 


C Kl,K2,Ml,Ll 


) 


CALL mind ( A ) 


100 

200 

300 


format 

format 

format 

RETURN 

ENO 


( *1* 5X *REOUCE* ) 
( 5X *IR* 15, «X 
( / 5X *STORAGE FOR 


All* IS, ax 
0* 17 / ) 


*12* IS,/ > 


SUBROUTINE INITAL ( Ml ,Ll,MD,M2,L0,L2,Nr)0F,N0EL. IB ) 
COMMON 0(33000) , JO (300) ,K0( 300) , MB (300) 

Mo^e establishes* 0-»LI,L2,ANO Jn,K0,MB— Ml,Mi 

LO a LI 


call lOBlN ( RMREAD, l,Ni)OF, I ) 
CALL lOBlN ( aHREAO,l,NOEL,l ) 
IB a NOEL /NQOF 


L2 s LI ♦ NOEL - ! 

LIMSL2.L1+1 

CALL IO0IN ( «HRtAO*l,l)(Ll).LlM ) 

M2 s MJ ♦ NDOF - i 
LI s LI - IB 

DO I MsmI,M2 

LI s LI ♦ IB 

JO(H) s LI 

KO(M) a H2 . IB 

I MB(M) a IB 

Ml a M2 ♦ 1 

LI a Ll « IB 

RETURN 
ENO 

SUBROUTINE CHOL« ( L» IS, ILS, ILT.mbaNO, IC ) 
COMMON D(33B00) , jO(30e) ,KD(3BB) ,MB(3»e) ,NUM 


DO a I=IL8 ,ILT,mbAN0 

L = L ♦ I 

M a IS^ 1 

ID a 1 ♦ IE 

Jl s I ♦ I 

J2 s I ♦ L 

D(I) a Om / O(I-ID) 

DO 2 JsJt,J2 

M a M ♦ 1 
IC = IS- KP(M) 

ID a I - JD(M)-IC 

Kl a I 

IF ( IC,LT,0 ) Kl a Kl-IC 

K2 a J .1 

SUM a 0,0 

DO 1 KSK1,K2 
C NIIM B NJM ♦ 1 

1 SUM 3 SUM ^ D(K) * D(K-IO) 

2 D(J) a ( D(J)-5UH ) / 0(1-10) 

SUM a 0,0 

DO 3 KaI,J2 

C NIJM B NUM ♦ 1 

3 SUM S SUM ♦ D(K)*D(K) 

<» 0(J2fl) a SORT ( 0(J2fl)-SUM ) 

return 

end 

SUBROUTINE UPDATE ( II, ! 2 , Ml, LI ) 

COMMON 0(33000) ,.10(300 ),K0( 300) , MB (300) 

13 s 11 1 

NUM a 0 
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DO 2 X3XWI2 


J 


XC ■ KO(X) - 13 

J1 * JO(I) 

If ( XCtLTa^ ) J1 • Jt « 1C 

J 2 » JO(l) ♦ M 8 (T) • 1 

DO 1 JsJWJ2 

NUM 3 NUN ♦ 1 
O(NUM) a D(J) 


t 8 X • 13 

JO(L) 8 NUM ♦ J1 • J2 

K0(L) 8 KO(X) - 13 

XF ( XC.LT.0 ) KO(L) s e 


XF ( IC.LT.B ) 


Ha(U 8 HB(X) 
M8CL) a H0 (l) 


X3 


2 CONTXNUC 


Ml 8 HI • 13 
ll 8 HUH ♦ 1 


RETURN 

ENO 

subroutine SOLVE ( NOPT»NOEQ#IO ) 

COMMON R(4323) 

CALL lOBXN ( «HREA0#2,R,N0E0 ) 

XF ( XO,GT,0 ) PRINT 100 
XF ( XO.GT.0 ) PRINT 200« (R(I),|al,NOEO) 
CALL second ( T 1 ) 

CALL FPASS ( NOPT»JE ) 


CALL TXMEX ( 0,0, T1 ) 

CALL SECOND ( T2 ) 

CALL BPA8S ( NOPT»JE ) 

CALL TIMEX ( 0»0,T2 ) 

XF ( XO,GT,0 ) PRINT 100 

IF ( XOfGT«0 ) PRINT 200, (R(I),Ial,NOEO) 

FRINT 300, T1 
PRINT aO0,T2 


100 

200 

300 

400 


FORMAT ( *1* 5X *SOLVE* ) 
format ( 5X *R * 20Ffl,t ) 

FORMAT ( *1* 5X *TlME FOR FPASS* F6.2 ) 
FORMAT ( 6X *TIME FOR 6PASS* F6,2 ) 


RETURN 

End 

subroutine FPaSS (N0PT,JE) 

COMMON R(4323) 

COMHON / TIUNS / 0( IM) ,NOFPE,N«OROS,H8*NO,LePJ( Ulll ) 

JE B tf 




00 3 N00F.S1 ,NQPT 


Call lOBlN ( 7HRf A0S«P»«,0» 1«3*LCPJ( NODE ) > 

P«!NT 200*NoFRE,NKO«DS#hH4NO 
20<« format ( 515 ) 

JE » JE ♦ NdFRF 

JS s JE • HBANO ♦ 1 

12 ■ CnoFRE • 1) * M8AN0 ♦ 1 

M a MBANO • NDFRE - 2 


00 3 XsU!2,HBAN0 

B « M ♦ I 
J2a I ♦ M 
IDs JS • I 
SUM B 0,0 


IF ( I.GT.J2 ) goto 2 
00 1 JaI,J2 

I SUM e SijH ♦ 0(J) * R(J«-ID) 


2 J a J2 ♦ I 

R(J+I0) a C R(jtlO)-SUM ) / 0(J) 


3 CONTINUE 

RETURN 

ENO 

SUproUTTnE brass fNOPT»JE) 

COMMON R(a323) 

COMMON / TRANS / O(ie 0 ),NOFRE,NwOROS#M 8 ANO,LCP.Uiaan 

00 NOOEal^NOPT 

LOC a NOPT-NOOEft 

CALL lOBIN ( 7HREA08«P,8,D,195,LCPJ ( LOC ) ) 

JT 8 JE 

JS a JT*PBANOfl 


00 2 lai^NOFRE 

Jl a (NOFRE • I) a mranO ♦ I 
J2 a J| ^ MbANO • 1 . 1 
P(JT) a R(JT) / 0(J2^n 
10 a JS - Jl 


IF (J1,CT,J2) GOTO 3 
DO 1 


IS7 


|g| 


0(J) * R(JT) 


1 R(J+IO) s R(J>ID) - 

2 JT s JT • I 


3 J£ s JE • NOFRE 

R CONTINUE 

RETURN 

ENO 

SUBROUTINE WIND ( NTAPE ) 

1 call IOBIN ( bHhRITER, NTAPE ) 

2 IP ( lOBlNiaHTEST, NTAPE)) 2,3,3 

3 CALL lOfliN ( 6HREWIN0, NTAPE ) 
RETURN 

ENO 
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APPENDIX C.3 

COLSOL COLUMN EQUATION SOLVER 


program test ( input, output, tape 5= I RPUT»TAPE6=0UTPUT,T ape t»UPE2, 

I TAPE'S, TAPEa,TAPE7,TAPE8) 

COMMUNJ S ,R(?H,5) , A (5d) ,MA (2H) ,KA C^J) ,0 (50) ,M9(2M) ,KB(a) , 

1 *120,5) ,MX(20) ,KX(a) ,KH(20) ,Ne(20) ,NEB( 10) ,NEBL( 10) , 

2 0(?P) ,n:T(6) 

MAXTsSO 

mAXC=20 
00 5 Ni=l 
5 ftiT (U ) SN 
NT (5) =7 
NT(6 )s0 

i« continue 

READ (5,ltiP»0) NEO,LEO,N0LK,NBLL,NLO 
IF (NFQ.EfJ.U) STOP 
HRiTEIb, IPMH) NEO,LEQ,NBLK,N0LL,NLO 
c read s upper triang full 

00 IPO 1=1, nEQ 

READ (5,2000) (S(J,1) ,J=1, n 
00 100 J=l,l 
100 S(I,J)=S(J,I) 

CALL PRMAT (S,NEO,NE5,20, IHS) 

DO 120 1=1, NLO 

120 READ (5,2000) (R(J,I),J=1,NE0) 

CALL PRMAT (R,NEQ,nlO,20, IHR) 

C PROFILE S 

00 200 1=1, NEC 
KH(I)=I 

DO 180 J=1,IM 
IF (S(I,J)) 200,170,200 
170 KH(I)=I«J 

180 continue 

200 CONTINUE 

*»RITE (6,1000) (kh(I),IsI,neq) 
read (5,1000) (NE8(I) , I=l,NBLK) 
wRTTE(6, 1 000) (MEB(1 ) , 1=1 ,N8LK) 

NL = 0 

CO 500 N=1,N0LK 
NFzNL+I 

nLsNF+ne8(N)-1 

KA(1 )=N 

KA(3)SNF 

KA(a)=.NL 

KA (2)=N 

LS0 


00 300 I=NF,NL 

NB(I)SN 

JF=I-KM(I)+1 

KHISNEO^JFfl 

00 250 I=JF,I 

L=L+l 

250 A(L)=S(J,n 

HA(I-NF+1)3L 

c find lo.»iest prev, block to operape 
IF (N.EQ.l) GO To i90 
KA(2)=MINM (KA (2) ,NH(NE0 -KhI+1) ) 

300 continue 

call tapes (NT(1),N,A,MAXT,MA,HAXC,KA,2) 
k^RlTE (6,3000) N,KA 

'^RITE (6,P0O0) (ma(I-nF^I),I = NF,NL) 

9000 FORMAT (UH ha ,1615) 
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BLANK 


o o o 


CALL PHAR (Af B»KA,KB»MA,PB) 

50B continue 

«EAD (5, lanq) (NEBL(n,Isi,N8LL) 
write (6, 1000) (NERL(I) , IsUnbLL) 

NL = 0 

00 6U0 N3l,NHLL 
NFsN|,tl 

NLaNF+NEBL(N)-l 

KACDsN 

KA(2)s0 

AA(3)snF 

KAC«)=NL 

Ls0 

uO 55 m IsNFjNL 
DO sail!) Jsj,neQ 
L=L^1 

5a0 A(L)=R(J,I) 

C 550 MAfl.NF+nsL 
550 MACI-NF4-US0 

CALL TAPES tNT( 3 ) A,MaXT»HA,maxC#KA,2) 
write (6,«000) n,KA 
CALL PRR (A,H4,kA,NEQ) 

600 CONTINUE 

IF (LEQ,EQ,NEO) GO TO 650 
C SET 50LN LEQ-t-l TO NEO ON TAPE nt(5) 

WRITE (6,9991) 

9991 format (30H SOLVE SUBSTRUCT DISPLACEMENTS ) 

KKKsl 

call S 0 RE(A,B,D,MArM 8 ,NE 0 ,NE 0 ,N 8 LK,N 3 LL»MAXT,MAXC»NT»KKK) 
C SET SUaSTR DISPLS ON TAPE NT( 5 ) 

00 630 Nal,N8LL 

CALL Tapes (NT(a),N,A,MAXT,MA,MAXC»KA,n 

NTAaKA(U)«KA(3)+l 

LEDOaLEQ+l 

ICaw 

La0 

00 620 l3l,NTA 
DO 610 KaLEQQ,NE(J 
LSL^I 

610 a(L)3A(K+IC) 

620 ICalC+NEQ 

CALL Tapes (NT(5) ,N,B,L,MBrMAXC»KA,2) 

630 CONTINUE 
650 CONTINUE 

WRITE (6,9992) 

9992 FORMAT (13H SOLVE SYSTEM ) 

KKKrl 
NTSaa 
KKKa3 
KKKsa 
KKKrS 
NTS83 
KKKB2 

CALL 80RE(A,8,0,MA,MB,NE0,LE0#NBLK,NBLL,MaXT,MAXC,NT,KKK) 

K K K a 5 
NTSr3 

IF (LEO, GE. NEO) GO TO 850 
IF (FKK,NE,5) GO TO 650 

SET SURSTR DISPLACEMENTS IN RED LOAD BLOCKS ON TAPE NT(9) 
DO 830 Nal,NBLL 


A*****«a* 


C HEAD HEO LOADS 

CALL tapes (NT(U),M,A,MAXT,HA,MAXC»KA,i) 

NTA=KA(a)-KA(3)fl 

LLsNTA*(NiEQ-LFQ) 

call Tapes (ntcs) ,m,h,u »^ip,maxCpK0mu 

ICS0 
JC = 0 

LEtJ«=LEQ+l 
DO IsUnTa 

DO 810 JsLEQQ,NE(3 
810 A(IC + Jl=8(JC->-J*LEQ) 

ICsIC+NEQ 
820 JCsJC+NEQ-LFG 

call Tapes (nt(6) a,maxt,ma,maxc,ka.2) 

830 CONTINUE 
NT(a)=NT(6) 

850 continue 

CALL SORE(A,B,Of MA,HB,NEQ,LEQ»N6Lt<»N8LLf MAXT»MAXC#NT,KKK) 
DO 600 N=l,NPUL 

CALL TAPES (NTS»N,A,MAXT,MA,MAXCrKA,l) 

WRITE (6»6000) N 

CALL PHR (a,ma,ka,NEQ) 

C FORM S^OISP 

ML=KA(a)-KA(3)+l 
DO 750 L=t»NL 
DO 7A0 1=1, NEO 
SSs0,0 

LCs(L-U*nE0 
00 730 J=1,NE0 
730 SS=SS+S(I» J)*A(J+LC) 

7tt0 X(I,L)=SS 
750 CONTINUE 

white (6»7000) N 

call PRHAT (X,NEQfNL»20r IHX) 

808 CONTINUE 
GO TO 10 

1000 format (1615) 

2000 FORMAT (16F5,0) 

3000 format (/6H block KA ,ai5l 

4000 format (/ilH LOAD 8L0CK ,15/40 KA ,415) 

6000 FORMAT (/15H SOLUTION BLOCK ,15) 

7000 format (/12H K*X,, .BLOCK ,15) 

END 

SUBROUTINE PRMAT ( A , NR , NC , MM , H ) 

DIMENSION A(MM,1) 

white (6,2000) H 
2000 format (/9H matrix ,A5) 

DO 200 nsi,nC, 0 
NLsNf7 

IF (NL.GT.NC) NLsNC 
DO 2id0 1 = 1, NR 

200 WHITE (6,100^0 I, (All , J) , JsN,NL) 

RETURN 

1000 FORMAT (I5,8FJ4,6) 

END 

SUBROUTINE PRR f.R,MR,KP,NEQ) 

DIMENSION R(l),MP(l),KH(4) 

NL = KR(4)-KR(3)-f 1 
11=0 

DO 100 N=1,NL 

WRITE ( 6 , 10 tnn N, (H(I+II),I=l,NFO) 

100 II=II+NEO 
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kETURN 

format (3H R ,I<i»,iaF9,2) 
end 

subroutine sore (a,b,o,ma,mb,neo»i.eo,nblk,nbll»maxt,maxc»nt,kkk) 

— --S0L0TI0N OR REDUCTION OF LINEAR EQUATIONS STORED OUT OF CORE IN 
compacted active column BLOCKS OF APPROXIMATELY MAXT LOCATIONS, 
programmed by E WILSON AND H OOVEY JAN 1976 
DIMENSION A(MAXT) »B(MAXT) ,0(NE0) ,MA (MAXC ) , MB (MAXC ) , NT ( <» ) 

DIMENSION KA(<0,KB(a) 

arrays a and 8 ARE WORKING STORAGE AREAS FOR BLOCKS OF COLUMNS 
OF THE COEFFICIENT MATRIX OR LOAD VECTORS, WHERE MAXC IS THE 
MAXIMUM number of columns OR VECTORS IN A BLOCK, MA AND MB ARE 
INTEGER arrays OF LOCATION OF DIAGONAL TERMS IN THE A OR B ARRAYS, 
THE 0 ARRAY STORES DIAGONAL TERMS OF REDUCED MATRIX , 

KA(n»KB(n- BLOCK NUMBER OF A OR B BLOCK 

KA(2),K(,2)- number of lowest BLOCK TO OPERATE ON THIS BLOCK, 
KA(3),KB(3)^ number OF FIRST COLUMN IN BLOCK 
KA(a),KB(tt)- number OF LAST COLUMN IN BLOCK 
NBLK number of BLOCKS OF COEFFICIENT MATRIX TERMS 

NBLL NUMBER OF BLOCKS OF LOAD VECTORS 

NEQ NUMBER OF EQUATIONS IN COEFFICIENT MATRIX 

LEQ NUMBER OF LAST EQUATION TO BE REDUCED IN COEFFICIENT MATRIX 

NT(1) TAPE NUMBEP for STORAGE OF BLOCKS OF THE COEFFICIENT MATRIX 

NT(3) TAPE NUMBER FOR STORAGE OF BLOCKS OF THE LOAD VECTORS 

NT(«) TAPE NUMBER FOR STORAGE OF BLOCKS OF THE DISPLACEMENT OR 

REDUCED LOAD VECTORS 

KKKsl complete SOLUTION • REQUIRES LEQsNEQ 

KKKs2 FOHioARO reduction of COEFFICIENT MATRIX AND LOAD VEC^'P * 

KKKs3 FORWARD REDUCTION OF COEFFICIENT MATRIX ONLY 

KKKsa FORWARD REDUCTION OF LOAD VECTORS ONLY 

KKKS5 BACKSURSTITUTION ONLY, IF LEQ IS LESS THAN NEQ SUBSTRUCTURE 
DISPLACEMENT MUST BE PREVIOUSLY CALCULATED 
GO TO (l00,lP0,10B»3P'B,3e«) KKK 
C— — -BLOCK-0V-0LOCK TRIANGULARIZATION OF MATRIX 
100 OU 2D0 Nsl,N8LK 

C I, move PREVIOUSLY REDUCED BLOCK TO B 

IF (N,EQ,n GO TO U0 

call BLKOP (A,B,D,MA,M3,KA,KB,NEQ,LEQ,MAXT»MAXC, 1) 

C 2, READ NEW block FROM TAPE NT(1) 

IIH CALL TAPES (NT( D »N, A,MAXT»MA,MAXC,KA, I ) 

C 3, OPERATE ON BLOCK N wjTH BLOCK M 
M s KAC2) 

IF (M,EQ.N) GO TO 160 
NM3N«1 

IF (M,EQ,NM) GO TO 190 
C READ BLOCK M from TAPE NTC2) 

120 call tapes (NT(2),M,0,MAXT,Mfl,MAXC,K0,l) 

190 CALL BLKOP (A,B,D,MA,H0,kA,K0,NEQ,LEQ,MAXT»MAXC»2) 

M 5 M + 1 

IF (M,LT,N) GO TO 120 
C <1, SELF-REDUCTION OF BLOCK N 

160 call BLKOP (A,A,0,MA,MA,KA,KA,NEQrLEQ,MAXT»MAXC»3) 

C 5t WRITE REDUCED BLOCK N ON TAPE NT(2) 

call Tapes fNT(2)rN»A»MAXT,MA,MAXC»KA,2) 

200 continue 

IF (KKK,E0,3) return 
C—— REDUCTION OF LOAD BLOCKS 
300 NTA » NT(3) 

NTS s NT(9) 

IF (KKK,NE,5) go to 310 
NTA a NT(9) 

NTB a NT(3) 
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c 

c 


c 


c 


DO Msl»NHtL 

1, READ LOAD Block m khom tape nta 
Call tapes (NTa,m,a,maxt,ma,maxc,Ka, l ) 

2. forward reduction of load block M 
IF (KKK,EQ,51 go to 33b 

DO S^,^ N51,NjBLK 

call TAPES(NT(^),NrB»MAxT,M0,N|AXC»KB»l) 

CALL BLKOP (A»B,n,HA,M8,KA,KB,NE';,LEQ,KAXT»0AXC»«) 
32« CONTINUE 


3, EACKWAPO REOIicTIOM OF LOAD Hi OCK M 
33« IF (KKK,£Q,/l) GO TO 3<?0 
IF (KKK,EQ,2) GO TO 39^ 

335 NsNBLK 

3ao CALL tapes (NT (2) , N, B» HAXT »M0,MAXC »K8, n 

CALL BLKQP (A,H,0,ma,mb,ka,KP,NEQ,LEQ,maxT»MAXC»5) 

N = N*1 


IF (N.GT.O) go to 59fc' 

«. WRITE RESULTS ON TAPE NTS 
39fc) CALL TAPES(NTP,m, a,maXT,Ha,MAXC»KA»2) 

CONTINUE 

return 

END 

subroutine tapes fNT»NR, X,MAXT,MC,HAXCtKFiKK) 

COMHON/TAPES/NTAPEUU) 

dimension X(maxT) ,mc (maxC) ,kf (<o 


SUBROUTINE TO READ OR WRITE BLOCK OF INFORMATION (X.MC.KF) 
WHICH IS RECORD NUMBER MR ON TAPE NT, 

LOR,LHK SHOULD BE SET PROPORTIONAL TO THE COST OF ONE 
DUMMY TAPE READ OR ONE TAPE BACKSPACE RESPECTIVELY 
NTAPE(I) CONTAINS THE CURRENT RECORD POSITION OF TAPE I 
nTAPECI) need not be initialized if THE TAPE WAS WRITTEN 
PY S/R tapes or if The FIRST RECORD TO HE READ IS RECORD 1 
LOR=l 


LBKsl 


IF (NR.NE.n GO TO 
5b Rewind NT 

NTAPE(nt) s 1 
9i3 LR = NR*LDR 


LKs(iv r APE(\n-NW) *lhk 

IF (LR.LT.LK) GO TO 50 
100 IF (ntaPE(NT)-nR) 

200 REAP (NT) 

NTAPE(NT) r NTAPE(NT)+1 
GO TO 100 
300 backspace NT 

NTAPE(NT) =NTAPE(NT)-t 
GO TO 100 

aO0 IF (KK,EO,n READ (NT) X,MC,KF 
IF (KK.EO.2) write (NT) X,MC,KF 
NTaPE(NT)sNTAPF (NT)+1 
RETURN 
End 

SUBROUTINE HLKOP ( A , h , D , m a , MB , K A , K B , NEO , L EQ . M A X T , M a XC » KK ) 
dimension a (1 ) ,b ( I ) ,ma M ) ,mh( I ) ,KA((i) ,KB(a) ,D( 1) 

GO TO ( lUO, 300, 3V0,5M0) ,KK 

c 

C 1. MOVE A ARRAY TO H 

C 


100 on IM lsl,tt 

no K0(I)=KA(I) 

00 120 Isi,HAXC 
120 MP(I)=ha(I) 
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00 13vl IsIjMAXT 

13(> BmsA(I) 

RETURN 

2, REDUCE BLOCK A BY BLOCK fl 

30P NTA = KA(R)*.KA(3Ut 
NTHsKB(a)-K0(3>fl 
KJs0 

DC 390 JsijNTA 
IF (KK.EQ.a) MA(J)sJ*NtQ 
JJSKA(3)+J«l 
IK (KK.EO.a) JJsNEO 
KHJ=MA(J)-KJ 
iLaNTB 

IF (KK,EQ,5) ILSJ 
Kls0 

DO 380 1=1, IL 
IJsKB(3)>I-l 
KMIsrt0(n-KI 
KHSMIN0 (KHI,KHJ-JJ4.II) 

KF=MB(I)-KHfl 
*<LsMfl(I)*MAX0 (l,n-LE(3) 

SS=0,0 

ij=MA(j),jj+n 
IF (KK,NE,3) go to 310 
IF (I.EO.J) GO TO 355 
310 IF (KF,GT,KL) go to 380 
IC=IJ-mh(I) 

DO 320 K=KF,KL 
320 SS=SS+9(K)*a(k+1C) 

A(IJ)SA(IJ)-SS 
GO TO 380 

355 ko=jj-khj+i-kf 

If (KF.GT.kl) go to 38B 
C REDUCE COLUMN BY ITSELF 

00 370 K=KF,kl 
AOsD(KOfK) 

IF (aD) 360,370,3^0 

360 TcA(K)/AO 

SSsSSfT*A(K) 

A(K)=T 

370 CONTINUE 

A(IJ)=A(IJ)-SS 
380 KI=M0(I) 

IF (KK.EO.O) 0(II)=B(Kn 
KJsMA(J) 

IF (KK.EO.a) CO TO 390 
0(JJ)=A(KJ) 

390 CONTINUE 

IF (KK.EQ,«,ANO,KB(«),EU,NEQ) GO TO 900 
RETURN 

C DIVIDE LOADS 0Y DIAGONALS 

400 KJS0 

DO 450 J=1,NTa 
mA(J)sJ*NEO 
KHJaMA(J)-KJ 
IFsNEO-KRJfl 

IF (IF.GT.LEO) GO TO 450 
DO 440 II=lF,LEO 
IF (DCID) 430,440,430 
93 > A(!l4KJ)=Af IUKJ)/D(in 


o o r> 


COMTINUE 
«5i) kJ = Ma(J) 

RETURN 

3, RACK-SUMSTITUTIOn 

Sfcij II=KH(y) 

NTA=K4(a)-KA(3)+l 
52K I5IT-K0(3)>J 

IE (I,LE,^> RETURN 
KFsl 

IF (I,GT,l) KF=My(I-l)4.l 
KL = H6(I)-MAXkKl ,XI-LEfJ) 

IF (KF,GT,KL) go to 5)80 

KHIsM0(I)«KF-M 

Kj30 

PO 55k) Jsl,NTA 

ICsKJ+(II-KhI)-KF+1 

AIJ=A(KJ+II) 

DO 540 KsKF,KL 

540 A(K+lC)=A(KfIC)-B(K)*AIJ 
550 KJ=J*NEQ 

580 ii=n-i 

GO TO 520 
ENO 

SUBROUTINE PRAB ( A , B , K A , KB , M A , MR ) 

DIMENSION A(1 ) ,Btn ,MAU ) fMB( 1) ,KA(4) ,KH(a) 
DIMENSION RR(1«) 
nTasKA (a)p.KA(3) + l 

KFsl 

00 400 l3l,NTA 

MAsKA(3)+I-1 

KL=MA(I) 

00 320 J=1,NA 
320 RRCJ)=O,0 
KHsKL-KF 
L=NA-KH 

00 350 KsKF,KL 
KKsK-KF 

350 RR(KK+L)=A(K) 

WRITE (8,2000) I, (RR(J) , Jsl ,NA) 

400 KFsKUfl 

RETURN 

1000 format (2H B,I3,14F9,2) 

2000 FORMAT (2H A,I3,14F9,2) 

ENO 
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